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.

Mathematica y la resolución numérica de problemas de valores iniciales

Mathematica incluye una orden de resolución numérica de problemas de valores iniciales. Ésta es:

NDSolve[{ ecuacion, valor inicial}, funcion incognita, {t, t0, t1 }]

Por ejemplo, para el problema de valores iniciales propuesto se obtiene:

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

Esta orden produce normalmente una salida como la anterior, en la que no aparece una solución explícita sino una función obtenida interpolando una serie de puntos que el método obtiene. Para realizar la gráfica de esta solución se ejecutan las siguientes órdenes:

solucion = x[t] /. sol

grafndsolve = 
      Plot[solucion // Evaluate, {t, 0, 1}, 
        PlotStyle -> {Dashing[{0.005, 0.02}]}];

Show[grafndsolve, grafexacta]

Ejercicios

1.- Aproxime  en 50 puntos igualmente espaciados del intervalo [1, 2], mediante el método de Euler,  la solución del p.v.i.
      
      		x' = - (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.- Resolver el problema de valores iniciales propuesto con la orden NDSolve y compárela gráficamente con las soluciones aproximadas halladas mediante el método de Euler y con la solución exacta.

5.- Resuelva el problema de la desintegración de 200 g. de átomos radioactivos durante el tiempo [0,1],
                                                  x' = -0.01 xx(0) = 200, 
      mediante:
      a) El método de Euler en 20 puntos igualemente espaciados.
      b) La orden exacta DSolve.
      c) La orden aproximada NDSolve.
      Compare las soluciones obtenidas gráficamente.

5.- Resuelva el problema de la desintegración de 200 g. de átomos radioactivos durante el tiempo [0,1],
                                                  x' = -3 xx(0) = 200, 
      mediante:
      a) El método de Euler en 20 puntos igualemente espaciados.
      b) La orden exacta DSolve.
      c) La orden aproximada NDSolve.
      Compare las soluciones obtenidas gráficamente.