6.1 EDO con Taylor de 3 términos con Python



1. Ecuaciones diferenciales ordinarias aproximadas con Taylor

Referencia: Rodríguez 9.1.1 ejemplo p335. Chapra 25.1.3 p731

EDO con Taylor 3términos gráfico animado

En los métodos con Taylor para Ecuaciones Diferenciales Ordinarias (EDO) se aproxima el resultado a n términos de la serie, para lo cual se ajusta la expresión del problema a cada derivada correspondiente.

La solución empieza usando la Serie de Taylor para tres términos ajustada a la variable del ejercicio:

y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h E = \frac{h^3}{3!} y'''(z) = O(h^3)
EDO Taylor 3 términos esquema gráfico

A partir de la expresión de y'(x) y el punto inicial conocido en x[i],se busca obtener el próximo valor en x[i+1] al avanzar un tamaño de paso h.

Se repite el proceso en el siguiente punto encontrado y se continua hasta alcanzar el intervalo objetivo del ejercicio.

En éstos métodos, la solución siempre es una tabla de puntos xi,yi que se pueden usar para interpolar y obtener una función polinómica, siendo el polinomio la solución aproximada de la EDO.

Se presenta un ejercicio para describir una solución paso a paso.



2. Ejercicio

Referencia: Rodríguez 9.1.1 ejemplo p335. Chapra 25.1.3 p731

Se requiere encontrar puntos de la solución en la ecuación diferencial usando los tres primeros términos de la serie de Taylor con h=0.1 y punto inicial x0=0, y0=1

\frac{dy}{dx}-y -x +x^2 -1 = 0

que con nomenclatura simplificada:

y'-y -x +x^2 -1 = 0

3. Desarrollo Analítico - paso a paso

Al despejar el valor de  y' de expresión del ejercicio,

y' = y -x^2 +x +1

se puede obtener y" al derivar una vez,

y'' = y' -2x + 1

para luego combinar las expresiones en

y'' = (y -x^2 +x +1) -2x + 1

simplificando:

y'' = y -x^2 -x +2

Ecuaciones que permiten estimar nuevos valores yi+1 para nuevos puntos  muestra distanciados en i*h desde el punto inicial siguiendo las siguientes expresiones de iteración:

y'_i = y_i -x_i^2 + x_i +1 y''_i = y_i -x_i^2 - x_i +2 y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h
EDO Taylor 3 términos, resultado en intervalo

Se empieza evaluando el nuevo punto a una distancia x1= x0+h del punto de origen con lo que se obtiene y1 , repitiendo el proceso para el siguiente punto en forma sucesiva.


itera = 0 , x0 = 0, y0 = 1

y'_0 = 1 -0^2 +0 +1 = 2 y''_0 = 1 -0^2 -0 +2 = 3 y_1 = y_{0} + h y'_0 + \frac{h^2}{2!} y''_0 y_1 = 1 + 0.1 (2) + \frac{0.1^2}{2!} 3 = 1.215 x_1 = 0 + 0.1
EDO  con Taylor de 3 términos itera=0 gráfica

itera = 1 , x = 0.1, y = 1.215

y'_1 = 1.215 - 0.1^2 + 0.1 +1 = 2.305 y''_1 = 1.215 - 0.1^2 - 0.1 +2 = 3.105 y_2 = 1.215 + 0.1 (2.305) + \frac{0.1^2}{2!} 3.105 = 1.461 x_2 = 0.1 + 0.1 = 0.2
EDO con Taylor 3 términos itera=1 gráfica

itera = 2 , x = 0.2, y = 1.461

y'_2 = 1.461 - 0.2^2 + 0.2 +1 = 2.621 y''_2 = 1.461 - 0.2^2 - 0.2 +2 = 3.221 y_3 = 1.461 + 0.1 (2.621) + \frac{0.1^2}{2!} 3.221 = 1.7392 x_3 = 0.2 + 0.1 = 0.3
EDO con Taylor 3 términos itera=2 gráfica

Completando los puntos con el algoritmo y realizando la gráfica se obtiene:

EDO con Taylor 3 términos
i,[xi, yi, d1yi, d2yi, término 1, término 2 ]
0 [0. 1. 0. 0. 0. 0.]
1 [0.1   1.215 2.    3.    0.2   0.015]
2 [0.2      1.461025 2.305    3.105    0.2305   0.015525]
3 [0.3        1.73923262 2.621025   3.221025   0.2621025  0.01610513]
4 [0.4        2.05090205 2.94923262 3.34923262 0.29492326 0.01674616]
5 [0.5        2.39744677 3.29090205 3.49090205 0.32909021 0.01745451]
>>>

Observación, note que los resultados de las derivadas, se encuentran desplazados una fila para cada iteración. Asunto a ser considerado en la gráfica de las derivadas en caso de incluirlas.

EDO Taylor 3 términos itera=muestras gráfico

Nota: Compare luego los pasos del algoritmo con el método de Runge-Kutta de 2do orden.



4. Algoritmo en Python para EDO con Taylor de 3 términos

Para simplificar los cálculos se crea una función edo_taylor3t() para encontrar  los valores para una cantidad de muestras distanciadas entre si h veces del punto inicial [x0,y0]

# EDO. Método de Taylor con3 términos 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
import numpy as np

# INGRESO
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1
# d1y = y', d2y = y''
d1y = lambda x,y: y - x**2 + x + 1
d2y = lambda x,y: y - x**2 - x + 2
x0 = 0
y0 = 1
h = 0.1
muestras = 5

# Algoritmo como función
def edo_taylor3t(d1y,d2y,x0,y0,h,muestras,
                 vertabla=False, precision=6):
    ''' solucion a EDO usando tres términos de Taylor,
    x0,y0 son valores iniciales, h es el tamaño de paso,
    muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+4),dtype=float)
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0,0,0]
    
    xi = x0 # valores iniciales
    yi = y0
    for i in range(1,tamano,1):
        d1yi = d1y(xi,yi)
        d2yi = d2y(xi,yi)
        yi = yi + h*d1yi + ((h**2)/2)*d2yi
        xi = xi + h
        
        term1 = h*d1yi
        term2 = ((h**2)/2)*d2yi
        
        tabla[i] = [xi,yi,d1yi,d2yi,term1,term2]
    if vertabla==True:
        np.set_printoptions(precision)
        print(' EDO con Taylor 3 términos')
        print('i,[xi, yi, d1yi, d2yi, término 1, término 2 ]')
        for i in range(0,tamano,1):
            print(i,tabla[i])
        
    return(tabla)
 
# PROCEDIMIENTO
tabla = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
n = len(tabla)
 
# SALIDA
print('EDO con Taylor 3 términos')
print('i,[xi, yi, d1yi, d2yi, término 1, término 2 ]')
for i in range(0,n,1):
    print(i,tabla[i])


5. Gráfica en Python para EDO con Taylor

Instrucciones adicionales para el algoritmo

# GRAFICA  --------------------
import matplotlib.pyplot as plt

titulo = 'EDO con Taylor 3 términos'
i = muestras # iteración en gráfica

titulo = titulo+', i='+str(i)
xi = tabla[:,0]
yi = tabla[:,1]
d1yi = tabla[:,2]
d2yi = tabla[:,3]

plt.plot(xi[0:i+2],yi[0:i+2]) # iteraciones
plt.plot(xi[0],yi[0],'o',
         color='red', label ='[x0,y0]')
plt.plot(xi[1:i+2],yi[1:i+2],'o',
         color='green', label ='[x[i],y[i]]')

if i<muestras: # gráfica para una iteración
    plt.plot(xi[i+1],yi[i+1],'o',color='orange',
             label ='[x[i+1],y[i+1]]')
    plt.plot(xi[i:i+3],yi[i:i+3],'.',color='gray')
    plt.plot(xi[i:i+2],[yi[i],yi[i]], color='orange',
             label='h',linestyle='dashed')
    term_1 = h*d1yi[i+1]
    plt.plot([xi[i+1],xi[i+1]],[yi[i],yi[i]+term_1],
             color='green',label='term2',linestyle='dashed')
    term_2 = (h**2)/2*d2yi[i+1]
    plt.plot([xi[i+1],xi[i+1]],[yi[i]+term_1,yi[i]+term_1+term_2],
             color='magenta',label='term3',linestyle='dashed')
# entorno de gráfica
plt.title(titulo)
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.tight_layout()
plt.show() # plt.show() #comentar para la siguiente gráfica

Tarea: Realizar el ejercicio con más puntos muestra, donde se visualice que el error aumenta al aumentar la distancia del punto inicial [x0,y0]



6. Cálculo de Error con la solución conocida

La ecuación diferencial ordinaria del ejercicio tiene una solución conocida, lo que permite encontrar el error real en cada punto respecto a la aproximación estimada.

y = e^x + x + x^2

Note que el error crece al distanciarse del punto inicial

Para las siguientes instrucciones, comente la última línea #plt.show() antes de continuar con:

# ERROR vs solución conocida
y_sol = lambda x: ((np.e)**x) + x + x**2

yi_psol = y_sol(xi)
errores = yi_psol - yi
errormax = np.max(np.abs(errores))

# SALIDA
print('Error máximo estimado: ',errormax)
print('entre puntos: ')
print(errores)

# GRAFICA [a,b+2*h]
a = x0
b = h*muestras+2*h
muestreo = 10*muestras+2
xis = np.linspace(a,b,muestreo)
yis = y_sol(xis)

plt.plot(xis,yis,linestyle='dashed', label='y solución conocida')
plt.legend()
plt.show()

Se puede observar los siguientes resultados:

Error máximo estimado:  0.0012745047595
entre puntos: 
[ 0.  0.000170  0.000377  0.000626  0.000922  0.00127 ]


Differential equations, a tourist's guide | DE1. 3Blue1Brown. 31-mayo-2019
Los Fundamentos de Ecuaciones Diferenciales que Nadie te Explica. 3Blue1Brown Español. 19-Junio-2024



Unidades MN