Métodos numéricos en ecuaciones diferenciales y MATLAB

Esta parte del blog está dedicada a los métodos numéricos aplicados a obtener aproximaciones de soluciones de ecuaciones diferenciales de primer orden. Para aclarar de modo más técnico este problema consideremos la ecuación diferencial con condición inicial

\frac{dy}{dx}=f(x,y), \hskip2cm y(a)=y_0

El problema que buscamos resolver consiste en determinar una aproximación para el valor y(b) para b\neq a.

 

Usted lector se preguntará seguramente, ¿porqué es necesario resolver este problema si ya conozco métodos que me permiten resolver muchas ecuaciones diferenciales?

Usted, amable lector de este blog, podrá imaginarse sin mucha dificultad ecuaciones diferenciales de primer orden, que al menos a primera vista  parezcan no tenener solución o al menos que hallar la solución no sea una tarea simple . Por ejemplo, pensemos en la ecuación diferencial.

 

\frac{dy}{dx}=(2x+y)x^2\sin^3(e^{x+y}).

Su solución pareciera no es fácil de determinar, en caso de que exista. Hablando de existencia de soluciones,   podemos responder afirmativamente a los problemas de existencia y unicidad de soluciones de la ecuaión diferencial para todo punto de en plano, usando el teorema de existencia y unicidad de ecuaciones diferenciales y que el lector encontrará en la sección correspondiente de este blog). Como es muy factible que no podamos resolver esta ecuación diferencia con condición inicial dada y(a)=y_0, es por ello que sería de gran utilidad hallar aproximaciones para valores y(b) para b\neq a.

 

En el documento siguiente, en formato pdf «SOLUCIÓN NUMÉRICA EDO’s», podrá hallar información sobre los métodos más elementales que se usan para aproximar soluciones para una ecuación diferencial de primer orden y para sistemas de ecuaciones diferenciales de primer orden. Desafortunadamente no contiene explicaciones explícitas del desarrollo de los distintos métodos, pero espero le sea de ayuda si ha tenido un primer contacto con este problema. Conforme el tiempo me lo vaya permitiendo, iré desarrollando en esta sección las ideas contenidas en el referido documento.

Contiene información sobre

* Método de Euler

* Método de Euler mejorado, también conocido como Predictor-Corrector.

* Método del punto medio

* Método de Runge-Kutta de cuarto orden

* Ideas de las estructuras computaciones que sirven de modelos para desarrollar programas.

* Programas en Matlab

* Ejemplos de los distintos métodos desarrollados

* Obtención del método de Runge-Kutta de segundo orden y determinación de algunos casos particulares de la familia de métodos de Runge-Kutta de segundo orden obtenidos.

* Compraciones numéricas de los métodos usados.

* Desarrollos sobre como extender los métodos de Euler y punto medio  para sistemas de ecuaciones diferenciales de primer orden.

 

Link para Documento

SOLUCION NUMERICA EDO’s

El siguiente link tiene un conjunto de problemas para resolver numéricamente usando los métodos que se han revisado

PROBLEMAS NUMÉRICOS

 

 

Estimado lector, permítame sugerirle la entrada

ECUACIONES DIFERENCIALES Y MATLAB

que se encuentra en este blog y que espero también le sea de ayuda.

APROXIMACIONES NUMÉRICAS DE SOLUCIONES DE ECUACIONES DIFERENCIALES MEDIANTE DESARROLLOS DE TAYLOR

 

Un tema que no se trata es el de la relación entre soluciones numéricas de una ecuación diferencial y el desarrollo de Taylor para una función.

Antes de ver dicha relación retomaremos aspectos que se usarán sobre el desarrollo de Taylor de una función y aproximaciones de Taylor de orden1,2 y 3 para una función dada.

APROXIMACIONES DE TAYLOR DE ORDENES 1,2 Y 3.

Recordemos que si en un punto x_0 conocemos información de todos los órdenes para una función y; es decir,  conocemos y(x_0), y^{(1)}(x_0), \dots,y^{(n)}(x_0),\dots entonces podemos obtener los valores de  y para puntos x^* suficientemente cercanos a x_0. Más precisamente,

y(x^*)=y(x_0)+y^{(1)}(x_0)\cdot (x^*-x_0)+\frac{1}{2!}y^{(2)}(x_0)\cdot (x^*-x_0)^2+\frac{1}{3!}y^{(3)}(x_0)\cdot (x^*-x_0)^3+\dots+\frac{1}{n!}y^{(n)}(x_0)\cdot (x^*-x_0)^n+\dots

Al tomar h=x^*-x_0, la expresión anterior toma la forma

y(x^*)=y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2+\frac{1}{3!}y^{(3)}(x_0)\cdot h^3+\dots\linebreak+\frac{1}{n!}y^{(n)}(x_0)\cdot h^n+\dots.

Otra forma de escribir la expresión anterior, cuya obtención queda fuera de los alcances de esta sección (el lector interesado puede consultar pràcticamente cualquier texto de Análisis Matemático) está dada por

y(x^*)=y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2+\frac{1}{3!}y^{(3)}(x_0)\cdot h^3+\dots\linebreak+\frac{1}{n!}y^{(n)}(x_0)\cdot h^n+ \frac{1}{n!}\int_{x_0}^{x^*}(x^*-t)^n y^{(n+1)}(t)dt.

donde el último sumando \frac{1}{n!}\int_{x_0}^{x^*}(x^*-t)^n y^{(n+1)}(t)dt se conoce como el residuo.Este residuo no es más que la diferencia entre y(x^*) y la expresión

y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2+\frac{1}{3!}y^{(3)}(x_0)\cdot h^3+\dots\linebreak+\frac{1}{n!}y^{(n)}(x_0)\cdot h^n.

Esta última expresión es conocida como la aproximación de orden n del valor y(x^*). Esto lo que escribimos como

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2+\frac{1}{3!}y^{(3)}(x_0)\cdot h^3\linebreak+\dots+\frac{1}{n!}y^{(n)}(x_0)\cdot h^n.

De esta forma, tenemos las aproximaciones de los órdenes siguientes

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h                                (aproximación de orden 1)

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2          (aproximación de orden 2)

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2+\frac{1}{3!}y^{(3)}(x_0)\cdot h^3          (aproximación de orden 3).

APROXIMACIONES NUMERICAS DE SOLUCIONES UANDO DESARROLLO DE TAYLOR DE ORDEN 1

A continuación mostraremos como dar la descripción del desarrollo de Taylor de orden 1

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h          (1)

para obtener una aproximación del valor y(x^*)  (con x^* suficientemente cercano a x_0) que toma la solución y(x) de la ecuación diferencial

y'=f(x,y)     (2)

que satisface la condiciónn inicial y(x_0)=y_0.

Observemos que de (1), para obtener la aproximación al valor  y(x^*) se conoce el valor y(x_0) es cual es igual a y_0, y solamente resta por conocer el valor y'(x_0).

Al ser y(x) solución de (2), satisface la igualdad

y'(x_0)=f(x_0,y(x_0))

pero y(x_0)=y_0. Así,

y'(x_0)=f(x_0,y(x_0))=f(x_0,y_0). Es decir, y'(x_0)=f(x_0,y_0), y al sustituir esta igualdad en (1), obtenemos

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h=y(x_0)+f(x_0,y_0)\cdot h,

con lo que se tiene la aproximación

y(x^*)\cong y(x_0)+f(x_0,y_0)\cdot h,,

la cual corresponde al método de Euler.

 

Es decir, el método de Euler es la aproximación para la solución de la ecuación diferencial que se obtiene usando el desarrollo de Taylor de orden 1.

 

APROXIMACIONES NUMERICAS DE SOLUCIONES UANDO DESARROLLO DE TAYLOR DE ORDEN $$

Del inciso anterior tenemos que al usar la igualdad  y'(x_0)=f(x_0,y_0) en el desarrollo de Taylor de orden 2

y(x^*)\cong y(x_0)+y^{(1)}(x_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2,

obtenemos

y(x^*)\cong y(x_0)+f(x_0,y_0)\cdot h+\frac{1}{2!}y^{(2)}(x_0)\cdot h^2,          (3)

y ahora solamente resta obtener una expresión para y^{(2)}(x_0), la cual se logra calcular a partir de la misma ecuación diferencial

y'(x)=f(x,y(x))

y el uso de la ragla de la cadena

Así,

y''(x_0)=\frac{df(x_0,y(x_0))}{dx}=\frac{\partial f \left(x_0,y(x_0)\right)}{\partial x}+\frac{\partial f\left(x_0,y(x_0)\right)}{\partial y} y'(x_0)=\frac{\partial f \left(x_0,y_0\right)}{\partial x}+\frac{\partial f\left(x_0,y_0\right)}{\partial y} f(x_0,y_0).

Es decir

y''(x_0)=\frac{\partial f \left(x_0,y_0\right)}{\partial x}+\frac{\partial f\left(x_0,y_0\right)}{\partial y} f(x_0,y_0).     (4).

 

Finalmente, sustituyendo (4) en (3), obtenemos

y(x^*)\cong y(x_0)+h\cdot f(x_0,y_0) +\frac{h^2}{2!}\left[\frac{\partial f \left(x_0,y_0\right)}{\partial x}+\frac{\partial f\left(x_0,y_0\right)}{\partial y} f(x_0,y_0)\right].       (5)

Esta expresión corresponde a la aproximación de la solución de la ecuación diferencial usando el desarrollo de Taylor de orden 2.

 

 

El algoritmo asociado a (5) para obtener la familia de aproximaciones \left\{y_i\right\}_{i=0,\dots, n} a los valores \left\{y(x_i)\right\}_{i=1,\dots, n} que la solución y(x) toma en los puntos a=x_0,x_1,\dots,x_n=b corresponde a

y_{i+1}= y(x_i)+h\cdot f(x_i,y_i) +\frac{h^2}{2!}\left[\frac{\partial f \left(x_i,y_i\right)}{\partial x}+\frac{\partial f\left(x_i,y_i\right)}{\partial y} f(x_i,y_i)\right].       (6)