Resolución numérica de problemas de valores iniciales

En esta práctica desarrollamos algunos métodos para construir soluciones numéricas  de problemas de valores iniciales asociados a una ecuación diferencial de primer orden. Concretamente, consideraremos el p.v.i. 

		x' = f(t, x),   x(a) = α, 
		
donde f es una función suficientemente regular definida en un dominio D y ((Subscript[t, 0]), (Subscript[x, 0])) ϵ D.  Supongamos que tiene una única solución definida en un intervalo I  que contiene al punto a,  y sea b  un valor tal que [a, b] ⊂ I.  Dividamos este intervalo en n ≥ 1 partes iguales mediante los nodos {(Subscript[t, j])}_(0 ≤ j ≤ n), donde

		(Subscript[t, j]) = a + j h,  j = 0, 1,  ... ,  n, 

siendo h = ((b - a)/n). Nos planteamos construir sucesiones {(Subscript[x, j])}_(0 ≤ j ≤ n) tales que (Subscript[x, j]) sea una aproximación del valor de la solución en el nodo (Subscript[t, j]), que es x((Subscript[t, j])).  Tales sucesiones dependerán de h.

El método de Euler

El método de Euler se escribe como

		(Subscript[x, j + 1]) = (Subscript[x, j]) + h f((Subscript[t, j]), (Subscript[x, j])),  j = 0, 
      1, ..., n - 1(Subscript[x, 0]) = α

Vamos a aplicarlo para resolver numéricamente el p.v.i. 

		x' = t + x,    x(0) = 1
		
en el intervalo [0, 1].

La programación del método de se lleva a cabo de forma que se puedan modificar todos los elementos que intervienen en el proceso.

Comenzamos subdividiendo el intervalo en 20 partes iguales. Con los nodos obtenidos se define la lista nodo y las aproximaciones de la solución del p.v.i. en aquéllos se almacenarán en la lista valor, que se inicializa con la condición inicial y se forma por etapas.

\(Clear[f, t, x, a, b, n, h, j, nodo, valor]
\)\(f[t_, x_] := t + x; 

a = 0.; 
b = 1.; 

α = 1.; 

    valor = {α}; 


\(n = 20; 
h = \(b - a\)\/n; 



    \(nodo = Table[a + j\ h, {j, 0, n}]; \)\)\)

Hallamos las aproximaciones y las incluimos sucesivamente en la lista valor.

\(\(For[j = 2, 
        j ≤ n + 1, \(j++\), 
	\ \ \ \ \ \ \ \ AppendTo[valor, 
          valor[[j - 1]] + 
            h\ f[nodo[[j - 1]], 
                valor[[j - 
                    1]]]]
\ \ \ \ \ \ \ ];\)\)

Representamos gráficamente los resultados obtenidos.

grafEuler1 = 
        ListPlot[
          Table[{nodo[[j]], 
              valor[[j]]}, {j, 
              n + 1}]];

La solución general de la ecuación diferencial del p.v.i se halla con el  comando DSolve.

En nuestro caso la orden precisa es la siguiente:

DSolve[{x'[t]==t + x[t],x[0] == 1},x[t],t]

Mostramos la gráfica de la misma.

grafexacta = Plot[-1 + 2 (E^t) - t, {t, a, b}];

Finalizamos mostrandos conjuntamente las gráficas de la solución exacta y de la aproximada.

Show[grafEuler1, grafexacta];

Se observa claramente el error inherente a la solución numérica. La disminución de h debería redundar en una mejor aproximación. Hagamos, pues, 40 partes iguales y veamos el resultado.

\(valor = {α}; 

n = 40; 

h = \(b - a\)\/n; 


    nodo = Table[a + j\ h, {j, 0, n}]; 


    For[j = 2, j ≤ n + 1, \(j++\), 
	\ \ \ \ \ \ \ \ 
      AppendTo[valor, 
        valor[[j - 1]] + 
          h\ f[nodo[[j - 1]], 
              valor[[j - 1]]]]

      \ \ \ \ \ \ \ ]; 


    grafEuler2 = 
      ListPlot[Table[{nodo[[j]], 
            valor[[j]]}, {j, n + 1}]]; 

    
Show[grafEuler2, grafexacta]; \)

El resultado es manifiestamente mejor.

Comparemos ambas soluciones numéricas.

Show[grafEuler1, grafEuler2];

Además de haber aproximado en más puntos, en cada uno de los nodos de índice par hay dos aproximaciones diferentes de la solución. La aproximación obtenida al dividir por dos el valor del paso es más precisa.

Métodos de Taylor

El método de Euler se establece a partir de la fórmula de Taylor de primer orden. Es razonable pensar en obtener una mejora si se considera un orden más elevado.

Por ejemplo, si el orden es dos, el método se escribe como

	(Subscript[x, j + 1]) = 
        (Subscript[x, j]) + h f((Subscript[t, j]), (Subscript[x, j])) + 
          (1/(2!)) (h^2)( ∂(f/∂)t ((Subscript[t, j]), 
                    (Subscript[x, j])) + ∂(f/∂)x ((Subscript[t, j]), 
                    (Subscript[x, j])) f((Subscript[t, j]), (Subscript[x, j]))),  j = 0, 1,  ... , 
      n - 1. 
		
Además, (Subscript[x, 0]) = α.

Es el método de Taylor de orden 2.

Vamos a aplicarlo al mismo p.v.i. de la sección anterior.

La lista de nodos sigue siendo nodo y la lista valor es la de las aproximaciones. Aunque tanto la ecuación diferencial como los valores de a, b y n son los mismos, volvemos a definirlos para hacer que esta sección sea independiente.

Clear[f, t, x, a, b, n, h, j, nodo, valor]\(

f[t_, x_] := t + x; 

a = 0.; 
b = 1.; 

α = 1.; 

    valor = {α}; 

n = 20; 
h = \(b - a\)\/n; 



    \(nodo = Table[a + j\ h, {j, 0, n}]; \)\)

Generamos las aproximaciones. Recuerde la notación para las derivadas parciales.

For[j = 2,j ≤ n + 1,j++,
	         aux=valor[[
                    j - 1]]+h f[
                      nodo[[j - 1]], 
                      valor[[j - 1]]]+
			        (1/2)(h^2)(Derivative[1, 0][f][nodo[[
                            j - 1]], 
                          valor[[j - 
                              1]]]+Derivative[0, 1][f][nodo[[
                              j - 1]], 
                            valor[[j - 
                                1]]] f[
                            nodo[[j - 
                                1]], 
                            valor[[j - 
                                1]]]);
	        AppendTo[valor, aux]
       ];

La lista valor contiene los valores calculados.

valor

Representamos los puntos correspondientes a los nodos considerados y las aproximaciones obtenidas.

grafTaylor = 
        ListPlot[
          Table[{nodo[[i]], 
              valor[[i]]}, {i, 
              n + 1}]];

Mostramos la gráfica de la solución exacta junto con la anterior. Observe que es inapreciable la diferencia entre ambas.

Show[grafTaylor, grafexacta];

Note cómo la exactitud de la fórmula de Taylor empleada hace que con 21 nodos igualmente espaciados se obtenga un resultado muy bueno (compare con los obtenidos con el método de Euler con 21 y 41 nodos también igualmente espaciados).

Métodos de Runge-Kutta

La utilización de los métodos de Taylor exige calcular algunas derivadas sucesivas de la solución del problema. Los métodos de Runge-Kutta surgen como alternativa. Con ellos se pretende mantener el orden del error de la fórmula de Taylor empleada para deducir el correspondiente método de Taylor pero sin calcular derivadas. No vamos a deducir ninguno de ellos, sino a dar directamente uno conocido como método de Runge-Kutta clásico.

Se parte de (Subscript[x, 0]) = α. Para cada j = 1, 2, ..., n se generan los valores

		\(TraditionalForm
      \`\(x\_\(j + 1\) = 
        x\_j + \(1\/6\) \(h(K\_1 + 2  K\_2 + 2  K\_3 + K\_4)\), \)\)
		
donde

		(Subscript[K, 1]) = f((Subscript[t, j]), (Subscript[x, j]))\(TraditionalForm
      \`K\_2 = f(t\_j + \(1\/2\) h, x\_j + \(1\/2\) h\ K\_1)\)\(TraditionalForm
      \`K\_3 = f(t\_j + \(1\/2\) h, x\_j + \(1\/2\) h\ K\_2)\)(Subscript[K, 4]) = f((Subscript[t, j]) + h, (Subscript[x, j]) + h (Subscript[K, 3]))

Como ejemplo, vamos a resolver numéricamente en [0, 3]  el p.v.i.
		x' = -t x + (t^3) (x^3)x(0) = (1/2)
dividiendo el intervalo en 100 partes iguales.
Actuamos como en los casos de Euler y Taylor.

\(Clear[f, t, x, a, b, n, h, j, nodo, valor, k1, k2, k3, k4]
\)\(f[t_, x_] := t\^3\ x\^3 - t\ x; 

a = 0.; 
b = 3.; 


    valor = {1\/2.}; 

n = 100; 

h = \(b - a\)\/n; 


    nodo = Table[a + j\ h, {j, 0, n}]; 


    For[j = 2, j ≤ n + 1, \(j++\), 
	\ \ \ \ \ \ \ 
      k1 = f[nodo[[j - 1]], 
          valor[[j - 1]]]; 

      	\ \ \ \ \ \ \ 
      k2 = f[nodo[[j - 1]] + h\/2, 
          valor[[j - 1]] + 
            \(h\ k1\)\/2]; 
	\ \ \ \ \ \ \ 
      k3 = f[nodo[[j - 1]] + h\/2, 
          valor[[j - 1]] + 
            \(h\ k2\)\/2]; 
	\ \ \ \ \ \ \ 
      k4 = f[nodo[[j - 1]] + h, 
          valor[[j - 1]] + h\ k3]; 

      	\ \ \ \ \ \ \ 
      AppendTo[valor, 
        valor[[j - 1]] + 
          1\/6\ h\ \((k1 + 2\ k2 + 2\ k3 + k4)\)]
\ \ \ \ \ \ \ ]; \)

Representamos los valores obtenidos.

grafRK4 = 
        ListPlot[
          Table[{nodo[[j]], 
              valor[[j]]}, {j, 
              n + 1}]];

Podemos resolver el problema de valores iniciales propuesto con la orden DSolve.

DSolve[{x'[t]+t x[t]==(t^3) x[t]^3,x[0] == (1/2)},x[t],t]

Representamos la solución exacta.

grafexacta = 
        Plot[(1/(Sqrt[1 + 3 (E^((t^2))) + (t^2)])), {t, a, b}];

Finalizamos mostrando la gráfica de las aproximaciones encontradas por aplicación del método de Runge-Kutta junto con la de la solución exacta

Show[grafRK4, grafexacta];

Los métodos de Euler y de Runge-Kutta clásico para sistemas

Vamos a desarrollar los método de Euler y Runge-Kutta clásico para el siguiente problema de valores iniciales:
		x' = t + x + 2  yy' = -t + sen x - yx(0) = 1,  y(0) = 0
		
Aproximamos su solución en el intervalo [0, 1], subdividiéndolo en 25 partes iguales.

Comenzamos definiendo la función correspondiente al sistema diferencial, los valores iniciales, los extremos del intervalo, el número de partes en que se subdivide, el paso del método y los nodos resultates.

\(f[{t_, x_, y_}] = {t + x + 2  y, \(-t\) + Sin[x] - y}; 


    valor = {{1., 0.}}; 

a = 0.; b = 1.; 

n = 20; 


    h = \(b - a\)\/n; 

nodo = Table[a + j\ h, {j, 0, n}]; \)

Generamos las aproximaciones numéricas y las almacenamos en la lista valor.

For[j = 2, 
        j ≤ n + 1, j++, 
AppendTo[valor, 
          valor[[j - 1]] + 
            h f[{nodo[[j - 1]], valor[[j - 1, 1]], 
                  valor[[j - 1, 2]]}]]];

valor // TableForm

ListPlot[valor];

Resolvamos ahora el problema de valores iniciales

	x' = x + 2  yy' = -t + sen x - y + zz' = y - (e^z)x(0) = 1,  y(0) = 0,  z(0) = 0

\(f[{t_, x_, y_, z_}] = {x + 2  y, \(-t\) + Sin[x] - y + z, y - Exp[z]}; 
    

valor = {{1., 0., 0.}}; 

a = 0.; b = 1.; 

n = 20; 


    h = \(b - a\)\/n; 

nodo = Table[a + j\ h, {j, 0, n}]; 


    For[j = 2, j ≤ n + 1, \(j++\), 
      
AppendTo[valor, 
        valor[\([j - 1]\)] + 
          h\ f[{nodo[\([j - 1]\)], valor[\([j - 1, 1]\)], 
                valor[\([j - 1, 2]\)], valor[\([j - 1, 3]\)]}]]]; 



    valor // TableForm\)

Como las aproximaciones numéricas corresponden a puntos de (ℝ^3), necesitamos una orden para poder efectuar la correspondiente representación. Para ello, cargamos un paquete específico de gráficos.

<< Graphics`Graphics3D`

Usamos la orden ScatterPlot3D, con algunas opciones (la lista completa de opciones se consigue con la orden Options[Scatterplot3D]).

ScatterPlot3D[valor, PlotJoined -> True, 
        PlotStyle -> Thickness[0.01],  
        AxesLabel -> {t, x, y}];

Vamos a resolver los dos problemas anteriores por medio del método de Runge-Kutta clásico de cuatro evaluaciones. Redefinimos las funciones para hacer esta parte independiente de la anterior y que se pueda modificar fácilmente.

\(RK[f_List, {t_, x_List}, α_List, {a_, b_}, n_] := 
      Module[{ord = Length[f], h = \(b - a\)\/n, valor = {α}, j, nodo, 
          k1, k2, k3, k4, asignaciones}, 
nodo = 
          Table[a + j\ h, {j, 0, n}]; 
For[j = 2, 
          j ≤ n + 1, \(j++\), 
asignaciones1 = 
            Flatten[{t -> nodo[\([j - 1]\)]} ⋃ 
                Table[x[\([i]\)] -> valor[\([j - 1, i]\)], {i, 
                    ord}]]; 
k1 = \((f /. 
                asignaciones1)\); 
asignaciones2 = 
            Flatten[{t -> nodo[\([j - 1]\)] + h\/2} ⋃ 
                Table[x[\([i]\)] -> 
                    valor[\([j - 1, i]\)] + \(1\/2\) h\ k1[\([i]\)], {i, 
                    ord}]]; 
k2 = \((f /. 
                asignaciones2)\); 
asignaciones3 = 
            Flatten[{t -> nodo[\([j - 1]\)] + h\/2} ⋃ 
                Table[x[\([i]\)] -> 
                    valor[\([j - 1, i]\)] + \(1\/2\) h\ k2[\([i]\)], {i, 
                    ord}]]; 
k3 = \((f /. 
                asignaciones3)\); 
asignaciones4 = 
            Flatten[{t -> nodo[\([j - 1]\)] + h} ⋃ 
                Table[x[\([i]\)] -> 
                    valor[\([j - 1, i]\)] + h\ k3[\([i]\)], {i, 
                    ord}]]; 
k4 = \((f /. 
                asignaciones4)\); 
AppendTo[valor, 
            valor[[j - 1]] + 
              1\/6\ h\ \((k1 + 2\ k2 + 2\ k3 + 
                    k4)\)]]; 
{nodo, 
          valor}
\ \ \ \ \ \ \ ]\)

\(f[{t_, x_, y_}] = {t + x + 2  y, \(-t\) + Sin[x] - y}; 


    valor = {{1., 0.}}; 

a = 0.; b = 1.; 

n = 20; 

h = \(b - a\)\/n; 
    

nodo = Table[a + j\ h, {j, 0, n}]; 


    For[j = 2, j ≤ n + 1, \(j++\), 
      
k1 = 
        f[{nodo[\([j - 1]\)], valor[\([j - 1, 1]\)], valor[\([j - 1, 2]\)]}]; 
      
k2 = 
        f[{nodo[\([j - 1]\)] + h\/2, 
            valor[\([j - 1, 1]\)] + \(1\/2\) h\ k1[\([1]\)], 
            valor[\([j - 1, 2]\)] + \(1\/2\) h\ k1[\([2]\)]}]; 
      
k3 = 
        f[{nodo[\([j - 1]\)] + h\/2, 
            valor[\([j - 1, 1]\)] + \(1\/2\) h\ k2[\([1]\)], 
            valor[\([j - 1, 2]\)] + \(1\/2\) h\ k2[\([2]\)]}]; 
      
k4 = 
        f[{nodo[\([j - 1]\)] + h, valor[\([j - 1, 1]\)] + h\ k3[\([1]\)], 
            valor[\([j - 1, 2]\)] + h\ k3[\([2]\)]}]; 
      
AppendTo[valor, 
        valor[[j - 1]] + 
          1\/6\ h\ \((k1 + 2\ k2 + 2\ k3 + k4)\)]]
\)\(valor // TableForm
\)ListPlot[valor];

Resolvamos ahora el problema de valores iniciales

		x' = x + 2  y
y' = -t + sen x - y + z
z' = y - (e^z)
x(0) = 1,  y(0) = 0,  z(0) = 0

\(f[{t_, x_, y_, z_}] = {x + 2  y, \(-t\) + Sin[x] - y + z, y - Exp[z]}; 
    

valor = {{1., 0., 1.}}; 

a = 0.; b = 1.; 

n = 20; 


    h = \(b - a\)\/n; 

nodo = Table[a + j\ h, {j, 0, n}]; 


    For[j = 2, j ≤ n + 1, \(j++\), 
      
k1 = 
        f[{nodo[\([j - 1]\)], valor[\([j - 1, 1]\)], valor[\([j - 1, 2]\)], 
            valor[\([j - 1, 3]\)]}]; 
      
k2 = 
        f[{nodo[\([j - 1]\)] + h\/2, 
            valor[\([j - 1, 1]\)] + \(1\/2\) h\ k1[\([1]\)], 
            valor[\([j - 1, 2]\)] + \(1\/2\) h\ k1[\([2]\)], 
            valor[\([j - 1, 3]\)] + \(1\/2\) h\ k1[\([3]\)]}]; 
      
k3 = 
        f[{nodo[\([j - 1]\)] + h\/2, 
            valor[\([j - 1, 1]\)] + \(1\/2\) h\ k2[\([1]\)], 
            valor[\([j - 1, 2]\)] + \(1\/2\) h\ k2[\([2]\)], 
            valor[\([j - 1, 3]\)] + \(1\/2\) h\ k2[\([3]\)]}]; 
      
k4 = 
        f[{nodo[\([j - 1]\)] + h, valor[\([j - 1, 1]\)] + h\ k3[\([1]\)], 
            valor[\([j - 1, 2]\)] + h\ k3[\([2]\)], 
            valor[\([j - 1, 3]\)] + h\ k3[\([3]\)]}]; 
      
AppendTo[valor, 
        valor[[j - 1]] + 
          1\/6\ h\ \((k1 + 2\ k2 + 2\ k3 + k4)\)]]
\)\(valor // TableForm
\)ScatterPlot3D[valor, PlotJoined -> True, 
      PlotStyle -> Thickness[0.01], AxesLabel -> {t, x, y}];

Ejercicios

1.- Aproxime  en 50 puntos igualmente espaciados del intervalo [1, 2], mediante los métodos de Euler, Taylor de orden 2 y Runge-Kutta clásico, la solución del p.v.i.
      
      		x' = 3 (t^2) + 2 - ((t x)/(1 + 39 (t^2)))
x(1) = (39/8)

2.- Utilice el comando Dsolve para resolver el problema anterior y dibuje en tres gráficas distintas la solución exacta y los puntos correspondientes a las aproximaciones obtenidas en el ejercicio 1.

3.- Disminuya el tamaño del paso en el ejemplo desarrollado al programar el método de Euler, obtenga nuevas aproximaciones y compárelas gráficamente con los valores de la solución exacta.

4.- Aplique el método de Runge-Kutta clásico con 25 puntos igualmente espaciados al ejemplo del ejercicio anterior, y aplique los métodos de Euler y de Taylor de orden 2 al ejemplo desarrollado al programar el de Runge-Kutta.