Análisis Numéricos.

1. Introducción

Al momento de aplicar las Matemática s a situaciones del mundo real nos encontramos a menudo con problemas que no pueden ser resueltos analíticamente o de manera exacta y cuya solución debe ser abordada con ayuda de algún procedimiento numérico.

El análisis numérico o cálculo numérico es la rama de las matemáticas encargada de diseñar algoritmos para simular aproximaciones de solución a problemas en análisis matemático. Se distingue del cómputo simbólico en que no manipula expresiones algebraicas, sino números. 

El objetivo principal del análisis numérico es encontrar soluciones «aproximadas» para problemas complejos.  

Aunque hay muchos tipos de métodos numéricos, todos comparten una característica común : Invariablemente los métodos numéricos llevan  a cabo un buen número de tediosos cálculos aritméticos, aquí entra en escena la computadora y los lenguajes de programación.

No siempre se encuentran programas hechos para la resolución de aproximaciones, de allí la importancia de que sepamos como podemos planear algunas soluciones.

• 1)Raices de ecuaciones.
• 2)Sistemas de ecuaciones algebráicas lineales.
• 3)Ajustes a Curvas
• 4)Ecuaciones Diferenciales Ordinarias.
  

Veamos algunos ejemplos.

2. Método de Newton Raphson

El método de Newton-Raphson, permite hallar una raíz de una ecuación no-lineal siempre y cuando se parte de una buena estimación inicial de la misma.

El esquema iterativo de Newton puede derivarse del desarrollo de Taylor de la función alrededor de la estimación inicial.

Supongamos que tenemos esta función:

Si la graficamos ( en https://www.geogebra.org/) vemos que la gráfica sería:

Bien, si ampliamos la zona de la intersección de los ejes , vemos:

Esto indica que existen unas raices para esta función, para x entre [0;5] (3 para ser precisos).

El método de Newton-Raphson es un método abierto, en el sentido de que su convergencia global no está garantizada. La única manera de alcanzar la convergencia es seleccionar un valor inicial lo suficientemente cercano a la raíz buscada. 

¿Como sabemos cuando estamos cerca de una raíz? 

Respuesta: Cuando el valor de f(x₁) cambia de signo respecto de f(x₂) indica que tenemos una raíz entre x₁ y x₂.

Así, se ha de comenzar la iteración con un valor razonablemente cercano al cero (denominado punto de arranque o valor supuesto). 

La relativa cercanía del punto inicial a la raíz depende mucho de la naturaleza de la propia función; si ésta presenta múltiples puntos de inflexión (como es nuestro caso) o pendientes grandes en el entorno de la raíz, entonces las probabilidades de que el algoritmo diverja aumentan, lo cual exige seleccionar un valor supuesto cercano a la raíz. Una vez se ha hecho esto, el método linealiza la función por la recta tangente en ese valor supuesto. La abscisa en el origen de dicha recta será, según el método, una mejor aproximación de la raíz que el valor anterior. Se realizarán sucesivas iteraciones hasta que el método haya convergido lo suficiente.

Veamos un código en C++ que intenta encotrar una aproximación a la Raiz de esta función.

2.1. Código Fuente

//Autor Desconocido
#include <iostream>
using namespace std;
//Seccion de #includes
#include <iostream>
#include <iomanip> // setprecision
#include <cmath>
//Seccion constantes simbolicas
#define PRECISION 10 //Cantidad maxima de decimales
#define MAX_ITERACIONES 100 //Cantidad maxima de iteraciones
#define INTERVALOS 6 //Intervalos posibles
//Prototipos
void tabula(double a, double b, int intervalos);	// Muestra un # tabulado de intervalos
double f(double x);	// Retorna el valor de la funcion evaluada en x
double f_derivada(double x); // Retorna la derivada de la funcion evaluada en x

int main()
{
    double a;
    double b;
    double tolerancia;	// Tolerancia
    double x0; // Primera aproximacion
    double x1; // Siguiente aproximacion
    double error;	// Diferencia entre dos aproximaciones sucesivas: x1 - x0
    int iteracion; // # de iteraciones
    bool converge = true;
    
    cout << setprecision(PRECISION);	// Se establece la precision
    cout << "\nCalculo de las raices de una funcion aplicando el metodo de Newton - Raphson\n";
    cout << "\nIngrese el intervalo inicial [a,b]:" << endl;
    
    // Se ingresa el intervalo
    cout << "\na = ";
    cin >> a;
    
    cout << "b = ";
    cin >> b;
    
    // Se tabulan los valores de f para INTERVALOS intervalos
    tabula(a, b, INTERVALOS);
    
    // Se pide elegir una aproximaci�n inicial
    cout << "\nEscoja el punto inicial adecuado:   x0 = ";
    cin >> x0;
    
    // Se pide ingresar la tolerancia
    cout << "Tolerancia = ";
    cin >> tolerancia;
    
    // Iteraciones
    
    // Se imprimen los valores de la primera aproximacion
    cout << "\nAproximacion inicial:\n";
    cout << "x0 = " << x0 << "\n"
        << "f(x0) = " << f(x0) << "\n"
        << "f'(x0) = " << f_derivada(x0) << endl;
    
    iteracion = 1;
    do {
        
        if (iteracion > MAX_ITERACIONES) {
            converge = false;	// Se sobrepas� la m�xima cantidad de iteraciones permitidas
            break;
            
        } else {
            x1 = x0 - f(x0) / f_derivada(x0); // C�lculo de la siguiente aproximaci�n
            error = fabs(x1 - x0);	// El error es la diferencia entre dos aproximaciones sucesivas
            
            // Se imprimen los valores de la siguiente aproximaci�n x1, f(x1), f_derivada(x1), error
            cout << "\nIteracion #" << iteracion << endl;
            cout << "x" << iteracion << " = " << x1 << "\n"
                << "f(x" << iteracion << ") = " << f(x1) << "\n"
                << "f'(x" << iteracion << ") = " << f_derivada(x1) << "\n"
                << "error = " << error << endl;
            
            // La diferencia entre dos aproximaciones sucesivas es tambi�n conocida como error.
            // La condici�n de terminaci�n consiste en que que el error debe ser <= que la tolerancia dada
            // Si se cumple la condici�n de terminaci�n, se ha encontrado la raiz aproximada buscada.
            if (error <= tolerancia) { // Condici�n de terminaci�n
                converge = true;
                break;
                
                // Si no se cumple el criterio de terminaci�n, se pasa a la siguiente iteraci�n
            } else {
                x0 = x1;
                iteracion++;
            }
        }
        
    } while (1);
    
    // Respuesta final
    if (converge) {
        cout << "\n\nPara una tolerancia de " << tolerancia << " la raiz de f es: " << x1 << endl;
        
    } else {
        cout << "\n\nSe sobrepas� la m�xima cantidad de iteraciones permitidas" << endl;
    }
    
    cin.get();
    cin.get();
    return 0;
}

void tabula(double a, double b, int intervalos)
{
    int puntos = intervalos + 1;
    
    double ancho = (b - a) / intervalos;
    
    cout << "\n\tx\t\tf(x) " << endl;
    for (int i = 0; i < puntos; i++) {
        cout << "\t" << a << "\t\t" << f(a) << endl;
        a = a + ancho;
    }
}
//Esta es la funci�n de nuestro ejemplo
double f(double x)
{
    return x * exp(cos(x)) / 1.5 - 1;
    //return exp(-x) + 3 * x - 3;
}
//Esta es la derivada de la funci�n de nuestro ejemplo
double f_derivada(double x)
{
    return exp(cos(x)) * (1 - x * sin(x)) / 1.5;
    //return -1 * exp(-x) + 3;
}

3. NUMÉRICA DE ECUACIONES NO LINEALES

Trataremos sobre uno de los problemas más vastos de la aproximación numérica la solución de ecuaciones no lineales analizado de diferentes maneras desde la óptica analítica y su interpretación geométrica.
En el campo de la tecnología principalmente en la ingeniería nos encontramos generalmente  con el siguiente problema  determinar las raíces de la ecuación f(x) = 0.
Como la teoría de la difracción de la luz se precisa de la siguiente función:

Su gráfica es:

Otra podría ser la expresión para determinar las orbitas planetarias se precisa de la ecuación llamada ecuación de Kepler:

x- a. sen(x) =b

con a y b tomando diversos valores.

Es decir f(x) puede ser una función de variable real x, como es un polinomio en x, o como una función trascendente es decir:

O una función trascendente  ( Una función trascendente es una función que no satisface una ecuación polinómica cuyos coeficientes sean a su vez polinomios)

Para dar solución a estos problemas existen distintos algoritmos o métodos para encontrar las raíces de f(x) = 0, pero debemos tener en cuenta que ninguno es general, pues   en otras palabras no existe un método que funcione con todas las ecuaciones perfectamente.
Pero sólo en un reducido caso será posible obtener  las raíces  exactas de f(x) = 0, es decir cuando se trata de f(x)  factorizable,  en tal sentido tenemos:

donde  r_i son las raices.

El método de punto fijo o de aproximaciones sucesivas es, junto con el de Bisección, uno de los primeros métodos que se utilizaron para resolver ecuaciones algebraicas y trascendentes. No obstante que en la actualidad existen otros métodos más eficientes, el de punto fijo se considera el más simple en sus principios y en él se pueden apreciar claramente todas las características de un método de aproximaciones sucesivas.

3.1. Bisección

¿Que es el método de bisección?
El método de bisección es un algoritmo de búsqueda de raíces que trabaja dividiendo el intervalo a la mitad y seleccionando el sub-intervalo que tiene la raíz.

Este es uno de los métodos más sencillos y de fácil intuición para resolver ecuaciones en una variable. 

Se basa en el teorema del valor intermedio, el cual establece que toda función continua f en un intervalo cerrado [a, b] toma todos los valores que se hallan entre  f (a) y f (b). 

Esto es que todo valor entre f (a) y f (b) es la imagen de al menos un valor en el intervalo [a, b]. En caso de que f (a) y f (b) tengan signos opuestos, el valor cero sería un valor intermedio entre f (a) y f (b), por lo que con certeza existe un p o punto en [a, b] que cumple f (p) = 0. De esta forma, se asegura la existencia de al menos una solución de la ecuación f (a) = 0.

3.2. Punto fijo

El método del punto fijo es un método abierto, también llamado de iteración de un punto o sustitución sucesiva, que reordena la ecuación de la forma en que x esté del lado izquierdo de la ecuación, para buscarla intersección entre la recta identidad y la curva g(x), como se muestran en los siguientes ejemplos.

Plantemos ahora la función identidad y=x y una función x=g(x) llamada iteradora ( luego veremos como se obtiene a partir de f(x). Si graficamos:

Observe que la raiz de f(x) se encuentra en el mismo valor de x donde ocurre la intersección entre la recta identidad en color verde y la función g(x) en color naranja. Se usa la linea vertical en color morado en x=raiz como referencia de lo indicado.

El método consiste en establecer un punto inicia x₀ para la búsqueda, que se usa para calcular el valor g(x₀).

En la siguiente iteración el nuevo valor para x es g(x₀), que se refleja en la recta identidad y nuevamente se usa para calcular g(x).

El resultado iterativo se muestra en la figura animada, donde se observa que el resultado es convergente.

Obtenido de: http://blog.espol.edu.ec/analisisnumerico/punto-fijo-concepto/

Vamos a aplicar el método en la siguiente función:

La gráfica de esta función es:

Si planteamos una ecuación :

Primer paso:

trabajamos matematicamente para dejar la expresión de la siguiente forma  X=g(x) , a g se la conoce como función iteradora.

Podemos tener varias alternativas....

a)

b)

c)

Si planteamos:

  y luego despejando x:

d)

vamos a usar la  (d) .

Segundo paso

¿Cual es la raiz? Vemos en la gráfica que es x₀=2 lo que produce y=0 (raiz)

Tercer Paso

G(x₀) = x₁

Cuarto Paso

El siguiente paso es comparar x₁ con x₀ , resultando dos alternativas:  Sí x₀ y x₁  son iguale, encotré la solución, fin del problema. Si son distintas continúo de la siguiente manera.

Debemos resaltar que la sucesión x₀, x₁, x₂, x₃;…x_k… se va acercando al valor de la raíz r₁, de manera que x_k se encuentra más cerca de r que x_k-1 o se van alejando de la raíz.

4. Integral de Riemann

La integra de una función definida se puede aproximar por el método de Riemann. Recordemos que la integral de una función es el área bajo la curva.

$\int_{x_{0}}^{x_{f}}f\left(x\right)dx\approx\sum_{i=0}^{i=n}f\left(x_i\right)\varDelta x$

cuanto más pequeño el delta más exacta será la aproximación.

Ejemplo

Para calcular la integral de la siguiente función polinómica

$f\left(x\right)=ax^{3}+bx^{2}+cx+d$

$\int_{x_{0}}^{x_{f}}f\left(x\right)dx\approx f\left(x_{0}\right)\Delta x+f\left(x_{1}\right)\Delta x+...+f\left(x_{8}\right)\Delta x$

Ejercicio

Para que un tren complete el recorrido entre dos localidades, la locomotora debe ejercer una fuerza variable con la posición, que responde a la siguiente función:

$f\left(x\right)=-0,0035x^2 +18x +10000$

En un gráfico fuerza en función de la posición, el área bajo la curva es el trabajo que es lo que queremos calcular o aproximar.

Se puede aproximar dicha área utilizando la integral de Riemann. Esto consiste en dividir el espacio de interés en intervalos iguales y suponer que cada uno es un rectángulo de altura igual al valor de la función al comienzo del segmento y de ancho igual a la longitud del intervalo (aproximación de orden 0). Cuanto menor sea la longitud del intervalo, mayor sera la exactitud de la aproximación.

Realizar un programa en C++ para determinar el trabajo realizado por la locomotora, para ir desde la localidad 1 hasta la localidad 2 a 5000 m de distancia, utilizando una aproximación de orden 0. Analizar la influencia de la longitud del intervalo.

Implementar una función que evalúe el valor de la fuerza para una posición dada y otra función que vaya sumando las sucesivas áreas.

4.1. ejemplo solución

/*
 *  Germán Andrés Xander 2021  *
 *
 *   Programa que calcula el trabajo realizado por una fuerza, de función conocida,
 *   aproximando el área bajo la curva mediante interpolación de orden 0. 
*/

#include <iostream>
using namespace std;

float aproximacion_orden_0 (float inicio, float fin, int pasos); //función que calcula el área bajo la curva
float fuerza_en_posicion (float posicion); // función que evalúa la función en la posición dada.

int main(){
    float trabajo;
    trabajo= aproximacion_orden_0(0,5000,500);
    
    cout<<"El trabajo realizado por la Fuerza desde los 0m hasta los 5000m es de "<<trabajo<<" Joules";
    
    return(0);
}

float aproximacion_orden_0 (float inicio, float fin, int pasos ){
    float incremento=(fin-inicio)/pasos;
    float area_orden_0 = 0;
    for (float posicion=inicio; posicion < fin; posicion += incremento){
        area_orden_0 += fuerza_en_posicion(posicion)*incremento;
    }
    return (area_orden_0);
}

float fuerza_en_posicion (float posicion){
    float fuerza;
    fuerza=-0.0035*posicion*posicion+18*posicion+10000; //Fuerza en función de la posición F= -0,0035x^2 + 18x + 10000
    return(fuerza);
}

5. Ejercicio ortodrómica

Objetivo

Escribir un programa en C++ que calcule la ortodrómica, dada la latitud y longitud de ambas ubicaciones.

Marco teórico

La ortodrómica (del griego orthos "recto" y dromos "carrera") es el camino más corto (distancia) entre dos puntos de la superficie terrestre; es el arco del círculo máximo que los une, menor de 180 grados. [Wikipedia]

Suponiendo un radio terrestre constante de 6371 km, la ortodrómica se puede calcular de la siguiente manera:

$d=r\Delta\hat{\sigma}$

donde $$$r$ es el radio de la Tierra y $\Delta\hat{\sigma}$ es el ángulo central entre ambos puntos, que se puede obtener utilizando la ley esférica del coseno:

$\Delta\hat{\sigma}=\arccos(\text{sen}\phi_{i}\:\text{sen}\phi_{f}+\cos\phi_{i}\:\cos\phi_{f}\:\cos\Delta\lambda)$

C++ cuenta con una librería que incluye funciones trigonométricas (cmath), pero por cuestiones didácticas, aproximaremos por series de Taylor.

$\text{sen}(x)=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n+1)!}x^{2n+1}=x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}-\cdots$

$\cos(x)=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n)!}x^{2n}=1-\frac{x^{2}}{2!}+\frac{x^{4}}{4!}-\cdots$

$\arccos(x)=\frac{\pi}{2}-\sum_{n=0}^{\infty}\frac{(2n)!}{4^{n}(n!)^{2}(2n+1)}x^{2n+1}$

Sugerencias

• Además de las funciones seno, coseno y arco coseno, utilice funciones para calcular potencia y factorial.
• Utilizando de 6 a 10 términos en las series, se obtiene una muy buena aproximación.
• Defina una variable con el valor de pi.
• Para mejorar la precisión, utilice variables del tipo double.
• Para apreciar el resultado, defina la cantidad de decimales que se muestran en pantalla con:
  

  cout.precision(15);

5.1. ejemplo solución coseno

/*
 *  Germán Andrés Xander 2017  *
 *
 *   Programa que calcula el coseno de un número utilizando la serie de Taylor 
*/

#include <iostream>
using namespace std;

int factorial(int j);
long double calculaTerm(long double x, int j);

#define PI 3.1415926535897932384

int main(void){
   long double deg;
   long double rad;
   int i;
   long double result = 1;

   cout<<"Ingrese x en grados: ";
   cin>>deg;

   /* convertir x a radianes */
   rad = deg * PI / 180;
   for(i = 1; i <= 10; ++i){
       int j = 2*i;
       long double term = calculaTerm(rad, j);
       // suma si i es par y resta si es impar, según la serie de taylor
        if(i%2==0){
            result=result+term;
        }else{
            result=result-term;
        }
   }
   cout.precision(15);
   cout<< "el coseno de "<<deg<<" es "<<result;
} 

int factorial(int j){
    int fact=1;
    for (int i=2; i<j+1;i++){
       fact=fact*i;
   }
   return fact;
}

long double calculaTerm(long double x, int j){
    long double taylor=1;
    for(int i=0;i<j;i++){
        taylor=taylor*x;
    }
    return taylor/factorial(j);
 }