Ecuaciones diferenciales, por medio de diferencias finitas, método de Euler y c++.

La aproximación para resolver ecuaciones diferenciales, por medio de diferencias finitas es un método antiguo utilizado para obtener la solución numérica de ecuaciones diferenciales, aplicación desarrollada por Euler en 1768 ( 1 )

Solución de problemas de ecuaciones diferenciales por el metodo de las diferencias finitas.

Cuando se tiene una ecuacion diferencial como la indicada,

dy/dx = y ,

con una condición de frontera: cuando x=1 es y =5 (se llama de frontera porque es el valor en un extremo).

La solución analitica se realiza como se indica, separando las variables,

dy/y = dx

y realizando la integracion,

dy/y = dx

ln y = x + c1

y = e^(x+c1) = e^x . e^c1

(el simbolo ^ , significa de exponencial, es el numero de Neper elevado a la x, siendo este ultimo el exponente)

y = c2 . e^x

de la condicion de frontera, la ecuacion esta sujeta a x=1 es y=5

aplicando a la solucion propuesta es 5=c2.e^1, c2=5/e

luego y = (5/e).e^x

 

 

Se quiere buscar un método numérico y aproximado para resolverla, en el siguiente dominio,

1 ≤ x ≤ 3

Vamos a subdividir el intervalo de x que va de x = 1 hasta x = 3 en intervalos iguales h, de tal manera que si tenemos n puntos (n puede ser tan grande como queramos), las subdivisiones son iguales y de acuerdo al siguiente grafico hay 6 puntos y 5 subdivisiones,

 

: : : : : : : : http://82.166.171.228:8080/publicaciones/euler/1.PNG

Nota : en el ultimo grafico los valores de 1 a 6, son las cantidades de puntos y no los valores de x, cuyo dominio esta definido como 1 ≤ x ≤ 3, es decir x1=1 y el ultimo 3.

N es el numero de puntos y (N-1) es el numero de subdivisiones,

h = (xn - x1)/ (N 1)

Si N = 11, que es numero arbitrario de puntos (podriamos haber adoptado 6 como en el grafico o 100)

1 ≤ x ≤ 3

es h = (3 -1)/ (11 1) = 0.2

es dy/dx = y

el concepto de derivada es, dy/dx = lim (yx+h yx) / h si h ---> 0

: : : : : : : http://82.166.171.228:8080/publicaciones/euler/2.PNG

La pendiente de la recta  es la tangente trigonometrica del angulo que forma la recta con respecto al eje de las x y esta es la derivada

dy / dx = lim (yi+1 yi) / h cuando h---> 0

y en forma aproximada, dy / dx = (yi+1 yi) / h

se propuso como problema dy/dx = y

luego (yi+1 yi) / h = yi

yi+1 yi = yi. h

yi+1 = yi. h + yi = yi.( h + 1) ,

este esquema de diferencias finitas, porque el intervalo h no tiende a cero, es decir es una aproximacion del limite.

Comencemos a generar soluciones con

x1= 1 es y1= 5 por la condición de frontera

Entonces la y2 = 5.(0.2 + 1)

La y3 = 5.(0.2 + 1). (0.2 + 1)

y4 = 5.(0.2 + 1). (0.2 + 1). (0.2 + 1) = 5.(0.2 + 1)^3

y5 = 5.(0.2 + 1)^4

y6 = 5.(0.2 + 1)^5 ... y11 = 5.(0.2 + 1)^10

 

Utilizacion del lenguaje c++, para obtener soluciones.

El siguiente programa esta escrito en lenguaje c++, con ++Dev-C, en el que se generan soluciones del anterior problema obteniendo el resultado en dos formas, un resultado se ve en la pantalla del ordenador y el otro en forma de un archivo de texto, que el lector puede consultar resultados, siendo yr, el valor exacto de la solucion.

El programa fue escrito con N= 51,

#include <cstdlib>

#include <fstream>

#include <iostream>

#include <stdio.h>

#include <math.h>

 

using namespace std;

int main(int argc, char *argv[])

 

{

 

// Este tipo de lineas de codigo que comienzan por '//' son comentarios

// esta es la ecuacion diferencial a resolver dy/dx = y , este es el valor de la subdivision h = (xn - x1)/ (N 1)

//Si N = 51 y el dominio es 1 <= x <= 3

// yi+1 = yi.( h + 1)

 

float h,x,y,xn, x1,y1,e ,yr;

int N, k ,sw;

e=2.71828 ;

sw=0;

N=51;

// el dominio es 1 <= x <= 3 , x es mayor e igual a 1 y x es menor e igual a 3

xn=3;

// valores de frontera para x1=1 es y1=5

x1 = 1;

y1 =5;

x = x1;

y = y1;

// valor de la subdivision

h = (xn-x1)/(N-1);

// Crea un fichero de salida

 

ofstream fs("resultados.txt");

cout << "la cantidad de puntos es N= "<< N << "-----valor de cada subdivision h="<< h << endl;

cout << " "<< endl;

cout << "x(1)="<< x1 << " y(1)="<< y1<<" h="<<h<< endl ;

 

cout << "----------------------------------------------------------------------" << endl;

fs << "la cantidad de puntos es N= "<< N << "-----valor de cada subdivision h="<< h << endl;

fs << "x(1)="<< x1 << " y(1)="<< y1<<" h="<<h << endl;

 

//las siguientes lineas son un metodo de resolucion de un problema,el algoritmo

 

for (k=2; k<=N ; k++)

{

if (sw=0) goto valor ;else goto continuar;

valor:

x = x1;

y = y1;

yr = (5/e)*pow(e,x);

cout << "x(1)="<< x1 << " y(1)="<< y1<< endl ;

cout << "x("<< 1 << ")=" << x<< " y(" << 1 <<")="<< y << " yr ="<< yr <<endl ;

sw=sw+1;

continuar:

y = y*(h+1) ; // ecuacion obtenida para el calculo

x=x+h;

// calculo del valor verdadero yr

yr = (5/e)*pow(e,x);

cout << "x("<< k << ")=" << x<< " y(" << k <<")="<< y << " yr ="<< yr <<endl ;

 

cout << "----------------------------------------------------------------------" << endl;

 

// Enviamos una cadena al fichero de salida:

fs << "x("<< k << ")=" << x<< " y(" << k <<")="<< y << " yr ="<< yr<< endl;

// Cerrar el fichero,

// para luego poder abrirlo para lectura:

}

fs.close();

system("PAUSE");

return EXIT_SUCCESS;

}

Referencias

[1] Los métodos numéricos.

Eduardo Ghershman, 28.4.2016