Algoritmos101

Cursos Python

Etiqueta: algoritmo Python

  • 7.2.1 EDP Elípticas método iterativo con Python

    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]

    ..

    1. Ejercicio

    Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

    Siguiendo el tema anterior en EDP Elípticas, para resolver la parte numérica asuma como:

    Valores de frontera: Ta = 60 , Tb = 25 , Tc = 50, Td = 70
    Longitud en:  x0 = 0, xn = 2 , y0 = 0, yn = 1.5
    Tamaño de paso dx = 0.25 dy = dx
    iteramax = 100 tolera = 0.0001

    Para la ecuación diferencial parcial Elíptica

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0


    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]
    ..

    2. EDP Elípticas: Método iterativo - Desarrollo Analítico

    A partir del resultado desarrollado en EDP elípticas:

    \frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda u_{i+1,j} +(-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +u_{i,j-1}=0

    y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

    se puede plantear que:

    u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

    que reordenando para un punto central desconocido se convierte a:

    u_{i,j} = \frac{1}{4} \big[ u_{i+1,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} \big]

    con lo que se interpreta que cada punto central es el resultado del promedio de los puntos alrededor del rombo formado en la malla.

    El cálculo numérico se puede realizar de forma iterativa haciendo varias pasadas en la matriz, promediando cada punto. Para revisar las iteraciones se controla la convergencia junto con un máximo de iteraciones.

    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]

    ..

    3. Algoritmo en Python

    Para el ejercicio dado, se presenta el resultado para tres iteraciones y el resultado final con un gráfico en 3D:

    u inicial:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
itera: 0   ;  u:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   53.12 51.41 50.98 50.87 50.84 50.84 44.27 25.  ]
 [60.   53.91 51.95 51.36 51.18 51.13 51.12 42.91 25.  ]
 [60.   54.1  52.14 51.5  51.3  51.23 51.21 42.59 25.  ]
 [60.   54.15 52.2  51.55 51.34 51.27 51.24 42.52 25.  ]
 [60.   58.85 58.07 57.72 57.58 57.52 57.5  48.76 25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
errado_u:  8.728094100952148
itera: 1   ;  u:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   53.83 51.69 50.98 50.75 50.68 49.02 41.73 25.  ]
 [60.   54.97 52.54 51.55 51.18 51.05 48.55 39.47 25.  ]
 [60.   55.31 52.89 51.82 51.39 51.23 48.4  38.85 25.  ]
 [60.   56.59 54.78 53.91 53.54 53.38 50.45 40.76 25.  ]
 [60.   61.17 60.91 60.6  60.42 60.33 57.38 48.29 25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
errado_u:  3.7443900108337402
itera: 2   ;  u:
[[50.   50.   50.   50.   50.   50.   50.   50.   50.  ]
 [60.   54.17 51.92 51.06 50.73 50.2  47.62 40.52 25.  ]
 [60.   55.5  52.97 51.76 51.23 50.3  46.45 37.7  25.  ]
 [60.   56.25 53.95 52.75 52.19 51.07 46.71 37.54 25.  ]
 [60.   58.05 56.71 55.9  55.47 54.33 49.8  40.16 25.  ]
 [60.   62.24 62.39 62.18 61.99 60.93 57.25 48.1  25.  ]
 [70.   70.   70.   70.   70.   70.   70.   70.   70.  ]]
errado_u:  2.0990657806396484
... continua
Método Iterativo EDP Elíptica
iteraciones: 41  converge =  True
xi: [0.   0.25 0.5  0.75 1.   1.25 1.5  1.75 2.  ]
yj: [0.   0.25 0.5  0.75 1.   1.25 1.5 ]
Tabla de resultados en malla EDP Elíptica
j, U[i,j]
6 [70. 70. 70. 70. 70. 70. 70. 70. 70.]
5 [60.   64.02 64.97 64.71 63.62 61.44 57.16 47.96 25.  ]
4 [60.   61.1  61.14 60.25 58.35 54.98 49.23 39.67 25.  ]
3 [60.   59.23 58.25 56.8  54.53 50.89 45.13 36.48 25.  ]
2 [60.   57.56 55.82 54.19 52.09 48.92 43.91 36.14 25.  ]
1 [60.   55.21 53.27 52.05 50.73 48.78 45.46 39.15 25.  ]
0 [50. 50. 50. 50. 50. 50. 50. 50. 50.]
>>>

    cuyos valores se interpretan mejor en una gráfica, en este caso 3D:

    EDP Elipticas Iterativo

    Instrucciones en Python

    # Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo
# d2u/dx2  + du/dt = f(x,y)
import numpy as np

# INGRESO
fxy = lambda x,y: 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson
# Condiciones iniciales en los bordes, valores de frontera
Ta = 60     # izquierda de la placa
Tb = 25     # derecha de la placa
#Tc = 50    # inferior 
fxc = lambda x: 50+0*x #  f(x) inicial inferior
# Td = 70   # superior
fxd = lambda x: 70+0*x #  f(x) inicial superior

# dimensiones de la placa
x0 = 0    # longitud en x
xn = 2
y0 = 0    # longitud en y
yn = 1.5
# discretiza, supone dx=dy
dx = 0.25  # Tamaño de paso
dy = dx

iteramax = 100  # maximo de iteraciones
tolera = 0.0001
verdigitos = 2   # decimales a mostrar en tabla de resultados
vertabla = True  # ver iteraciones

# coeficientes de U[x,t]. factores P,Q,R
lamb = (dy**2)/(dx**2)
P = lamb       # izquierda P*U[i-1,j]
Q = -2*lamb-2  # centro Q*U[i,j]
R = lamb       # derecha R*U[i+1,j]

# PROCEDIMIENTO
# iteraciones en xi,yj ancho,profundidad
xi = np.arange(x0,xn+dx/2,dx)
yj = np.arange(y0,yn+dy/2,dy)
n = len(xi)
m = len(yj)

# Matriz u[xi,yj], tabla de resultados
u = np.zeros(shape=(n,m),dtype=float)

# llena u con valores en fronteras
u[0,:]   = Ta   # izquierda
u[n-1,:] = Tb   # derecha
fic = fxc(xi)   # Tc inferior
u[:,0]   = fic
fid = fxd(xi)   # Td superior
u[:,m-1] = fid  
# estado inicial dentro de la placa
# promedio = (Ta+Tb+Tc+Td)/4
promedio = (Ta+Tb+np.max(fic)+np.max(fid))*(1/-Q)
u[1:n-1,1:m-1] = promedio

np.set_printoptions(precision=verdigitos)
if vertabla == True:
    print('u inicial:')
    print(np.transpose(u))

# iterar puntos interiores
U_3D = [np.copy(u)]
U_error = ['--']
itera = 0
converge = False
while (itera<=iteramax and converge==False):
    itera = itera +1
    Uantes = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            u[i,j] = P*u[i-1,j]+R*u[i+1,j]+u[i,j-1]+u[i,j+1]
            u[i,j] = u[i,j] - fxy(xi[i],yj[j])*(dy**2)
            u[i,j] = (1/-Q)*u[i,j]
    diferencia = u - Uantes
    errado_u = np.max(np.abs(diferencia))
    if (errado_u<tolera):
        converge=True
    U_3D.append(np.copy(u))
    U_error.append(np.around(errado_u,verdigitos))
    
    if (vertabla==True) and (itera<=3):
        np.set_printoptions(precision=verdigitos)
        print('itera:',itera-1,'  ; ','u:')
        print(np.transpose(u))
        print('errado_u: ', errado_u)
    if (vertabla==True) and (itera==4):
        print('... continua')
        
# SALIDA
print('Método Iterativo EDP Elíptica')
print('iteraciones:',itera,' converge = ', converge)
print('errado_u: ', errado_u)
print('xi:', xi)
print('yj:', yj)
print('Tabla de resultados en malla EDP Elíptica')
print('j, U[i,j]')
for j in range(m-1,-1,-1):
    print(j,u[:,j])

    La gráfica de resultados requiere ajuste de ejes, pues el índice de filas es el eje x, y las columnas es el eje y. La matriz y sus datos en la gráfica se obtiene como la transpuesta de u

    # GRAFICA 3D ------
import matplotlib.pyplot as plt
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x

fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U,
                       color ='blue')
graf_3D.plot(Xi[1,0],Yi[1,0],U[1,0],'o',color ='orange')
graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='red',
             label='U[i,j]')
graf_3D.plot(Xi[1,2],Yi[1,2],U[1,2],'o',color ='salmon')
graf_3D.plot(Xi[0,1],Yi[0,1],U[0,1],'o',color ='green')
graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')

graf_3D.set_title('EDP elíptica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('y')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

    EDP Elípticas [ concepto ] [ ejercicio ]
    Método iterativo: [ Analítico ] [ Algoritmo ]

     

  • 7.2 EDP Elípticas

    EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]
    ..

    1. EDP Elípticas

    Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

    Las Ecuaciones Diferenciales Parciales tipo elípticas semejantes a la mostrada:

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

    (ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)

    se interpreta como una placa metálica de dimensiones Lx y Ly, delgada con aislante que recubren las caras de la placa, y sometidas a condiciones en las fronteras:

    Lx = dimensión x de placa metálica
Ly = dimensión y de placa metálica
u[0,y]  = Ta
u[Lx,y] = Tb
u[x,0]  = Tc
u[x,Ly] = Td

    Placa Metalica 01

    Para el planteamiento se usa una malla en la que cada nodo corresponden a los valores u[xi,yj]. Para simplificar la nomenclatura se usan los subíndices i para el eje de las x,  j para el eje t, quedando u[i,j].

    EDP Elipticas Iterativo

    vista superior, plano xy:

    Placa Metalica 02

    La ecuación se cambia a la forma discreta, usando diferencias divididas que se sustituyen en la ecuación, por ejemplo:

    \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^2}=0

    Se agrupan los términos de los diferenciales:

    \frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda u_{i+1,j} - 2\lambda u_{i,j} + \lambda u_{i-1,j} + u_{i,j+1} - 2 u_{i,j} + u_{i,j-1} = 0 \lambda u_{i+1,j} + (-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +u_{i,j-1}=0

    Obteniendo así la solución numérica conceptual. La forma de resolver el problema determina el nombre del método a seguir.

    Una forma de simplificar la expresión, es hacer λ=1, es decir  \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1, se determina que los tamaño de paso Δx, Δy son iguales.

    u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

    EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]

  • 7.1.2 EDP Parabólicas método implícito con Python

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método implícito: [ Analítico ] [ Algoritmo ]

    ..

    1. Ejercicio

    Referencia:  Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodríguez 10.2.4 p415

    Siguiendo el tema anterior en EDP Parabólicas, se tiene que:

    \frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}

    Ecuación Diferencial Parcial Parabólica método Implicito animado

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método implícito: [ Analítico ] [ Algoritmo ]

    ..

    2. Desarrollo Analítico

    En éste caso se usan diferencias finitas centradas y hacia atrás; la línea de referencia es t1:

    \frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \frac{\partial u}{\partial t} = \frac{u_{i,j} - u_{i,j-1} }{\Delta t}

    La selección de las diferencias divididas corresponden a los puntos que se quieren usar para el cálculo, se observan mejor en la gráfica de la malla.

    Barra Metalica 04

    Luego se sustituyen en la ecuación del problema, obteniendo:

    \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = K\frac{u_{i,j}-u_{i,j-1}}{\Delta t}

    De la gráfica se destaca que en la fórmula, dentro del triángulo solo hay DOS valores desconocidos, destacados por los punto en amarillo.
    En la ecuación se representa por U[i,j] y U[i+1,j]. Por lo que será necesario crear un sistema de ecuaciones sobre toda la línea de tiempo t1 para resolver el problema.

    EDP M Implicito 02

    Despejando la ecuación, se agrupan términos constantes: λ = \frac{\Delta t}{K (\Delta x)^2} .

    \lambda u_{i-1,j} + (-1-2\lambda) u_{i,j} + \lambda u_{i+1,j} = -u_{i,j-1}

    Los parámetro P, Q y R se determinan de forma semejante al método explícito:

    P = \lambda Q = -1-2\lambda R = \lambda Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j} = -u_{i,j-1}

    Los valores en los extremos son conocidos, para los puntos intermedios  se crea un sistema de ecuaciones para luego usar la forma Ax=B y resolver los valores para cada u(xi,tj).

    Por ejemplo con cuatro tramos entre extremos se tiene que:
    indice de tiempo es 1 e índice de x es 1.

    i=1,j=1

    Pu_{0,1} + Qu_{1,1} + Ru_{2,1} = -u_{1,0}

    i=2,j=1

    Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -u_{2,0}

    i=3,j=1

    Pu_{2,1} + Qu_{3,1} + Ru_{4,1} = -u_{3,0}

    agrupando ecuaciones y sustituyendo valores conocidos:
    \begin{cases} Qu_{1,1} + Ru_{2,1} + 0 &= -T_{0}-PT_{A}\\Pu_{1,1} + Qu_{2,1} + Ru_{3,1} &= -T_{0}\\0+Pu_{2,1}+Qu_{3,1}&=-T_{0}-RT_{B}\end{cases}

    que genera la matriz a resolver:

    \begin{bmatrix} Q && R && 0 && -T_{0}-PT_{A}\\P && Q && R && -T_{0}\\0 && P && Q && -T_{0}-RT_{B}\end{bmatrix}

    Use alguno de los métodos de la unidad 3 para resolver el sistema y obtener los valores correspondientes.

    Por la extensión de la solución es conveniente usar un algoritmo y convertir los pasos o partes pertinentes a funciones.

    Tarea: Revisar y comparar con el método explícito.

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método implícito: [ Analítico ] [ Algoritmo ]

    ..

    Algoritmos en Python

    Para la solución con el método implícito, se obtiene el mismo resultado en la gráfica y tabla. Aunque el algoritmo es diferente.

    EDP Parabolica
    algunos valores:

    iteración j: 1
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-40. -25. -25. -25. -25. -25. -25. -25. -35.]
resultado en j: 1 
 [31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]
iteración j: 2
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-46.01 -26.03 -25.18 -25.03 -25.01 -25.01 -25.08 -25.44 -37.57]
resultado en j: 2 
 [35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39]
iteración j: 3
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-50.25 -27.49 -25.55 -25.12 -25.03 -25.05 -25.24 -26.07 -39.39]
resultado en j: 3 
 [38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72]
EDP Parabólica - Método implícito 
lambda:  0.25
x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 10 ], primeras iteraciones
10 [60.   47.43 37.58 31.3  27.98 26.64 26.63 27.83 30.43 34.63 40.  ]
9 [60.   46.75 36.68 30.56 27.48 26.3  26.33 27.47 30.04 34.33 40.  ]
8 [60.   45.96 35.69 29.8  27.01 26.   26.06 27.12 29.6  33.99 40.  ]
7 [60.   45.01 34.6  29.02 26.57 25.73 25.81 26.76 29.13 33.58 40.  ]
6 [60.   43.87 33.39 28.24 26.16 25.51 25.58 26.41 28.6  33.09 40.  ]
5 [60.   42.45 32.06 27.48 25.8  25.32 25.39 26.07 28.03 32.48 40.  ]
4 [60.   40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40.  ]
3 [60.   38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40.  ]
2 [60.   35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40.  ]
1 [60.   31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>>

    Instrucciones en Python

    # EDP parabólicas d2u/dx2  = K du/dt
# método implícito
# Referencia: Chapra 30.3 p.895 pdf.917
#       Rodriguez 10.2.5 p.417
import numpy as np

# INGRESO
# Valores de frontera
Ta = 60
Tb = 40
#T0 = 25 # estado inicial de barra
fx = lambda x: 25 + 0*x # función inicial para T0
a = 0  # longitud en x
b = 1

K = 4     # Constante K
dx = 0.1  # Tamaño de paso
dt = dx/10

n = 100 # iteraciones en tiempo
verdigitos = 2   # decimales a mostrar en tabla de resultados

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx/2,dx)
fi = fx(xi)
m  = len(xi)
ultimox = m-1
ultimot = n-1
# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j = 0
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = -1 -2*lamb
R = lamb

# Calcula U para cada tiempo + dt
tj = np.arange(0,(n+1)*dt,dt)
j = 1
while not(j>=n):
    # Matriz de ecuaciones
    k = m-2
    A = np.zeros(shape=(k,k), dtype = float)
    B = np.zeros(k, dtype = float)
    for f in range(0,k,1):
        if (f>0):
            A[f,f-1]=P
        A[f,f] = Q
        if (f<(k-1)):
            A[f,f+1]=R
        B[f] = -u[f+1,j-1]
    B[0] = B[0]-P*u[0,j]
    B[k-1] = B[k-1]-R*u[m-1,j]
    # Resuelve sistema de ecuaciones
    C = np.linalg.solve(A, B)

    # copia resultados a u[i,j]
    for f in range(0,k,1):
        u[f+1,j] = C[f]

    # siguiente iteración
    j = j + 1

    # muestra 3 primeras iteraciones
    np.set_printoptions(precision=verdigitos)
    if j<(3+2): 
        print('iteración j:',j-1)
        print('A:\n',A)
        print('B:\n',B)
        print('resultado en j:',j-1,'\n',C)
        
# SALIDA
print('EDP Parabólica - Método implícito ')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:len(xi)],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',ultimox,'], primeras iteraciones')
for j in range(ultimox,-1,-1):
    print(j,u[:,j])

    Para realizar la gráfica se aplica lo mismo que en el método explícito

    # GRAFICA ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
plt.show()

    Sin embargo para la gráfica en 3D se ajustan los puntos de referencia agregados por las diferencias finitas.

    # GRAFICA en 3D
# vector de tiempo, simplificando en tramos
tramos = 10
salto = int(n/tramos)
if (salto == 0):
    salto = 1
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)

# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(tramos,m),dtype=float)
for k in range(0,tramos,1):
    U[k,:] = u[:,k*salto]
    tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)

fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U,
                       color ='blue',label='EDP Parabólica')
graf_3D.plot(Xi[2,0],Yi[2,0],U[2,0],'o',color ='orange')
graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')
graf_3D.plot(Xi[2,2],Yi[2,2],U[2,2],'o',color ='salmon')
graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='green')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

    Queda por revisar la convergencia y estabilidad de la solución a partir de los O(h) de cada aproximación usada. Revisar los criterios en Chapra 30.2.1 p891, Burden 9Ed 12.2 p727, Rodríguez 10.2.2 p409 .

    Tarea o proyecto: Realizar la comparación de tiempos de ejecución entre los métodos explícitos e implícitos. La parte de animación funciona igual en ambos métodos.  Los tiempos de ejecución se determinan usando las instrucciones descritas en el enlace: Tiempos de Ejecución en Python

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método implícito: [ Analítico ] [ Algoritmo ]

  • 7.1.1 EDP Parabólicas método explícito con Python

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método explícito: [ Analítico ] [ Algoritmo ]
    ..

    1. Ejercicio

    Referencia:  Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodríguez 10.2 p406

    Siguiendo el tema anterior en EDP Parabólicas, para resolver la parte numérica considere

    Valores de frontera: Ta = 60, Tb = 40, T0 = 25
    Longitud en x a = 0, b = 1,  Constante K= 4
    Tamaño de paso dx = 0.1, dt = dx/10

    Para la ecuación diferencial parcial parabólica:

    \frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}

    Ecuación diferencial parcial Parabolica animado ..

    2. Método Explícito

    Se usan las diferencias divididas, donde se requieren dos valores para la derivada de orden uno y tres valores para la derivada de orden dos:

    \frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \frac{\partial u}{\partial t} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t}

    La selección de las diferencias divididas corresponden a los puntos que se quieren usar para el cálculo, se observan mejor en la gráfica de la malla. La línea de referencia es el tiempo t0

    Barra Metalica 03

    Luego se sustituyen en la ecuación del problema, obteniendo:

    \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = K\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

    De la gráfica se destaca que en la fórmula solo hay un valor desconocido, destacado por el punto en amarillo dentro del triángulo. En la ecuación se representa por U[i,j+1].

    Despejando la ecuación, se agrupan términos constantes:

    λ = \frac{\Delta t}{K (\Delta x)^2}

    quedando la ecuación, con los términos ordenados de izquierda a derecha como en la gráfica:

    u_{i,j+1} = \lambda u_{i-1,j} +(1-2\lambda)u_{i,j}+\lambda u_{i+1,j}

    Al resolver se encuentra que que cada valor en un punto amarillo se calcula como una suma ponderada de valores conocidos, por lo que el desarrollo se conoce como el método explícito.

    La ponderación está determinada por los términos P, Q, y R.

    \lambda = \frac{\Delta t}{K(\Delta x)^2} P = \lambda Q = 1-2\lambda R = \lambda u_{i ,j+1} = Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j}

    Fórmulas que se desarrollan usando un algoritmo y considerando que al disminuir los valores de Δx y Δt la cantidad de operaciones aumenta.

    EDP Parabolica 01

    Queda por revisar la convergencia y estabilidad de la solución a partir de los O(h) de cada aproximación usada.

    Revisar los criterios en:  Chapra 30.2.1 p891 pdf915, Burden 9Ed 12.2 p727, Rodríguez 10.2.2 pdf409 .

    \lambda \leq \frac{1}{2}

    Cuando λ ≤ 1/2 se tiene como resultado una solución donde los errores no crecen, sino que oscilan.
    Haciendo λ ≤ 1/4 asegura que la solución no oscilará.
    También se sabe que con λ= 1/6 se tiende a minimizar los errores por truncamiento (véase Carnahan y cols., 1969).

    Para resolver la parte numérica considere que:

    Valores de frontera: Ta = 60, Tb = 40, Tc = 25 
Longitud en x: a = 0, b = 1 
Constante K = 4 
Tamaño de paso dx = 0.1, dt = dx/10

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método explícito: [ Analítico ] [ Algoritmo ]
    ..

    3. Algoritmo en Python

    Se muestra una tabla resumida de resultados a forma de ejemplo.

    Método explícito EDP Parabólica
lambda:  0.25
x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 5 ], primeras iteraciones de  100
5 [60.   44.21 32.93 27.29 25.41 25.05 25.18 25.98 28.4  33.23 40.  ]
4 [60.   42.77 31.29 26.37 25.14 25.   25.06 25.59 27.7  32.62 40.  ]
3 [60.   40.86 29.38 25.55 25.   25.   25.   25.23 26.88 31.8  40.  ]
2 [60.   38.12 27.19 25.   25.   25.   25.   25.   25.94 30.62 40.  ]
1 [60.   33.75 25.   25.   25.   25.   25.   25.   25.   28.75 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>>

    Se presenta la propuesta del algoritmo EDP parabólicas para el método explícito.

    # EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
import numpy as np

# INGRESO
# Valores de frontera
Ta = 60  # izquierda de barra
Tb = 40  # derecha de barra
#Tc = 25 # estado inicial de barra en x
fxc = lambda x: 25 + 0*x # f(x) en tiempo inicial

# dimensiones de la barra
a = 0  # longitud en x
b = 1
t0 = 0 # tiempo inicial, aumenta con dt en n iteraciones

K = 4     # Constante K
dx = 0.1  # muestreo en x, tamaño de paso
dt = dx/10

n = 100  # iteraciones en tiempo
verdigitos = 2   # decimales a mostrar en tabla de resultados

# coeficientes de U[x,t]. factores P,Q,R, 
lamb = dt/(K*dx**2)
P = lamb       # izquierda P*U[i-1,j]
Q = 1 - 2*lamb # centro Q*U[i,j]
R = lamb       # derecha R*U[i+1,j]

# PROCEDIMIENTO
# iteraciones en x, longitud
xi = np.arange(a,b+dx/2,dx)
m = len(xi) 
ultimox = m-1
ultimot = n-1

# u[xi,tj], tabla de resultados
u = np.zeros(shape=(m,n+1), dtype=float)

# u[i,j], valores iniciales
u[0,:] = Ta  # Izquierda
u[ultimox,:] = Tb  # derecha
# estado inicial de barra en x, Tc
fic = fxc(xi) # f(x) en tiempo inicial
u[1:ultimox,0] = fic[1:ultimox]

# Calcula U para cada tiempo + dt
tj = np.arange(t0,(n+1)*dt,dt)
for j  in range(0,n,1):
    for i in range(1,ultimox,1):
        # ecuacion discreta, entre [1,ultimox]
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]

# SALIDA
j_mostrar = 5
np.set_printoptions(precision=verdigitos)
print('Método explícito EDP Parabólica')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:j_mostrar],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',j_mostrar,'], primeras iteraciones de ',n)
for j in range(j_mostrar,-1,-1):
    print(j,u[:,j])

    Si la cantidad de puntos aumenta al disminuir Δx y Δt, es mejor disminuir la cantidad de curvas a usar en el gráfico para evitar superponerlas. Se usa el parámetro 'salto' para indicar cada cuántas curvas calculadas se incorporan en la gráfica.

    # GRAFICA ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
plt.show()

    Note que en la gráfica se toma como coordenadas x vs t, mientras que para la solución de la malla en la matriz las se usan filas y columnas. En la matriz el primer índice es fila y el segundo índice es columna.

    En el caso de una gráfica que se realice usando los ejes x,t,U en tres dimensiones se requiere usar las instrucciones:

    # GRAFICA en 3D ------
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)

# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(m,tramos),dtype=float)
for k in range(0,tramos,1):
    U[:,k] = u[:,k*salto]
    tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)
U = np.transpose(U)

fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U, color ='blue')
graf_3D.plot(xi[0],tk[1],U[1,0],'o',color ='orange')
graf_3D.plot(xi[1],tk[1],U[1,1],'o',color ='green')
graf_3D.plot(xi[2],tk[1],U[1,2],'o',color ='green')
graf_3D.plot(xi[1],tk[2],U[2,1],'o',color ='salmon',
             label='U[i,j+1]')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método explícito: [ Analítico ] [ Algoritmo ]

    Solving the heat equation | DE3. 3Blue1Brown. 16-junio-2019

    EDP Parabólica [ concepto ] [ ejercicio ]
    Método explícito: [ Analítico ] [ Algoritmo ]

  • 7.1 EDP Parabólicas

    EDP Parabólica [ concepto ] Método [ explícito ] [ implícito ]
    ..

    1. EDP Parabólicas

    Referencia:  Chapra 30.2 p.888 pdf.912, Burden 9Ed p714, Rodríguez 10.2 p406

    Las Ecuaciones Diferenciales Parciales tipo parabólicas, semejantes a la mostrada, representa la ecuación de calor para una barra aislada sometida a fuentes de calor en cada extremo.

    La temperatura se representa en el ejercicio como u[x,t]

    \frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}

    Barra Metalica 01

    EDP Parabolica Ani 01

    o con vista en 3D:

    EDP Parabolica 02

    Para la solución numérica, cambia la ecuación a su forma discreta usando diferencias finitas divididas que se sustituyen en la ecuación,

    \frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \frac{\partial u}{\partial t} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t}

    con lo que la ecuación continua se convierte a discreta:

    \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = K\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

    Para interpretar mejor el resultado, se usa una malla que en cada nodo representa la temperatura como los valores u[xi,tj].

    Para simplificar nomenclatura se usan los subíndices i para el eje de las x y j para el eje t, quedando u[ i , j ].

    Barra Metalica 02

    En el enunciado del problema habían establecido los valores en las fronteras:

    • temperaturas en los extremos Ta, Tb
    • la temperatura inicial de la barra T0,
    • El parámetro para la barra K.

    El resultado obtenido se interpreta como los valores de temperatura a lo largo de la barra luego de transcurrido un largo tiempo. Las temperaturas en los extremos de la barra varían entre Ta y Tb a lo largo del tiempo.

    EDP Parabolica 01

    Tomando como referencia la malla, existirían algunas formas de plantear la solución, dependiendo de la diferencia finita usada: centrada, hacia adelante, hacia atrás.

    EDP Parabólica [ concepto ] Método [ explícito ] [ implícito ]

    But what is a partial differential equation? | DE2

    3Blue1Brown. 2019 Abril 21

    Tarea: Revisar ecuación para difusión de gases, segunda ley de Fick.ley de Fick

    La difusión molecular desde un punto de vista microscópico y macroscópico.

  • Unidad 7 EDP - Ecuaciones Diferenciales Parciales

    Parabólicas: [ explícito ] [ implícito ]

    Elípticas: [ iterativo ] [ implícito ]

    Hiperbólicas

    EDP Parabolica Ani 01

    Ecuación diferencial parcial Parabolica animado

  • 6.5 Métodos EDO con gráficos animados en Python

    animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

    Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para la solución de Ecuaciones Diferenciales Ordinarias, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.

    animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

    ..

    EDO con Taylor de 3 términos

    Edo con Taylor de 3 términos GIF animado

    Tabla de resultados:

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

    Instrucciones en Python

    # 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

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,6),dtype=float)
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        d1yi = d1y(x,y)
        d2yi = d2y(x,y)
        y = y + h*d1yi + ((h**2)/2)*d2yi
        x = x + h
        
        term1 = h*d1yi
        term2 = ((h**2)/2)*d2yi
        
        tabla[i] = [x,y,d1yi,d2yi,term1,term2]
    if vertabla==True:
        titulo = ' [xi,     yi,     d1yi,   d2yi,'
        titulo = titulo + '   término 1,   término 2 ]'
        np.set_printoptions(precision)
        print(' EDO con Taylor 3 términos')
        print(titulo)
        print(tabla)
        
    return(tabla)

# PROGRAMA -----------------
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1
# INGRESO.
# 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

# PROCEDIMIENTO
tabla = edo_taylor3t(d1y,d2y,x0,y0,h,muestras,
                     vertabla=True)

# SALIDA
##print(' EDO con Taylor 3 términos')
##print(' [xi,     yi,     d1yi,',
##      '   d2yi,   término 1,   término 2 ]')
##print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
xi = tabla[:,0]
yi = tabla[:,1]
plt.plot(xi,yi)
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
#plt.show() #comentar para la siguiente gráfica


# GRAFICA CON ANIMACION --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation

unmetodo = 'Edo con Taylor 3 términos'
narchivo = 'EdoTaylor3t' # nombre archivo GIF
muestras = 51 

# Puntos para la gráfica
a = xi[0]
b = xi[1]
term1 = tabla[:,4]
term2 = tabla[:,5]
dfi = tabla[:,2]
n = len(xi)

# Parametros de trama/foto
retardo = 1000   # milisegundos entre tramas
tramas = len(xi)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
ymax = np.max([yi[0],yi[2]])
ymin = np.min([yi[0],yi[2]])
deltax = np.abs(xi[2]-xi[0])
deltay = np.abs(yi[2]-yi[0])
graf_ani.set_xlim([xi[0],xi[2]+0.05*deltax])
graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.05*deltay])
# Lineas y puntos base
linea_fx, = graf_ani.plot(xi, yi,color='blue',
                          linestyle='dashed')
puntof, = graf_ani.plot(xi[0], yi[0],'o',
                        color='green',
                        label='xi,yi')
puntoa, = graf_ani.plot(xi[0], yi[0],'o',
                        color='Blue')
puntob, = graf_ani.plot(xi[1], yi[1],'o',
                        color='orange')

linea_h, = graf_ani.plot(xi, xi, color='orange',
                         label='h',
                         linestyle='dashed')
linea_term1, = graf_ani.plot(xi, xi,
                             color='green',label="h*y'[i]",
                             linestyle='dashed')
linea_term2, = graf_ani.plot(xi, yi, linewidth=4,
                             color='magenta',
                             label="((h**2)/2!)*y''[i]")
# Aproximacion con tangente
b0 = yi[0] - dfi[1]*xi[0]
tangentei = dfi[1]*xi + b0
linea_tang, = graf_ani.plot(xi, tangentei, color='dodgerblue',
                             label="tangente",
                             linestyle='dotted')

# Cuadros de texto en gráfico
txt_i  = graf_ani.text(xi[0], yi[0]+0.03*deltay,'[x[i],y[i]]',
                       horizontalalignment='center')
txt_i1 = graf_ani.text(xi[1], xi[1]+0.03*deltay,'[x[i+1],y[i+1]]',
                       horizontalalignment='center')
# Configura gráfica
graf_ani.axhline(0, color='black')  # Linea horizontal en cero
graf_ani.set_title(unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Nueva Trama
def unatrama(i,xi,yi,dfi,term1,term2):
    
    if i>1:
        ymax = np.max([yi[0:i+2]])
        ymin = np.min([yi[0:i+2]])
        deltax = np.abs(xi[i+1]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    else:
        ymax = np.max([yi[0],yi[2]])
        ymin = np.min([yi[0],yi[2]])
        deltax = np.abs(xi[2]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    # actualiza cada punto
    puntoa.set_xdata([xi[i]]) 
    puntoa.set_ydata([yi[i]])
    puntob.set_xdata([xi[i+1]]) 
    puntob.set_ydata([yi[i+1]])
    puntof.set_xdata([xi[0:i]]) 
    puntof.set_ydata([yi[0:i]])
    # actualiza cada linea
    linea_fx.set_xdata(xi[0:i+1])
    linea_fx.set_ydata(yi[0:i+1])
    linea_h.set_xdata([xi[i],xi[i+1]])
    linea_h.set_ydata([yi[i],yi[i]])
    linea_term1.set_xdata([xi[i+1],xi[i+1]])
    linea_term1.set_ydata([yi[i],yi[i]+term1[i+1]])
    linea_term2.set_xdata([xi[i+1],xi[i+1]])
    linea_term2.set_ydata([yi[i]+term1[i+1],
                           yi[i]+term1[i+1]+term2[i+1]])
    
    b0 = yi[i] - dfi[i+1]*xi[i]
    tangentei = dfi[i+1]*xi + b0
    linea_tang.set_ydata(tangentei)
    # actualiza texto
    txt_i.set_position([xi[i], yi[i]+0.03*deltay])
    txt_i1.set_position([xi[i+1], yi[i+1]+0.03*deltay])

    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_tang,
            linea_term1,linea_term2,
            txt_i,txt_i1,)
# Limpia Trama anterior
def limpiatrama(): 
    puntoa.set_ydata(np.ma.array(xi, mask=True))
    puntob.set_ydata(np.ma.array(xi, mask=True))
    puntof.set_ydata(np.ma.array(xi, mask=True))
    linea_h.set_ydata(np.ma.array(xi, mask=True))
    linea_term1.set_ydata(np.ma.array(xi, mask=True))
    linea_term2.set_ydata(np.ma.array(xi, mask=True))
    linea_tang.set_ydata(np.ma.array(xi, mask=True))
    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_tang,
            linea_term1,linea_term2,
            txt_i,txt_i1,)

# Trama contador
i = np.arange(0,tramas-1,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (xi,yi,dfi,term1,term2),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

    animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

    ..

    Runge Kutta de 2do Orden para primera derivada

    EDO Runge-Kutta 2do orden primera derivada _animado

     EDO con Runge-Kutta 2do Orden primera derivada
 [xi,     yi,     K1,    K2 ]
[[0.       1.       0.       0.      ]
 [0.1      1.2145   0.2      0.229   ]
 [0.2      1.459973 0.23045  0.260495]
 [0.3      1.73757  0.261997 0.293197]
 [0.4      2.048564 0.294757 0.327233]
 [0.5      2.394364 0.328856 0.362742]]

    Instrucciones en Python

    # EDO. Método de Runge-Kutta 2do Orden primera derivada 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras, vertabla=False, precision=6):
    ''' solucion a EDO con Runge-Kutta 2do Orden primera derivada,
        x0,y0 son valores iniciales
        muestras es la cantidad de puntos a calcular con tamaño de paso h.
    '''
    # Runge Kutta de 2do orden
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+2),dtype=float)
    
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (1/2)*(K1+K2)
        xi = xi + h
        
        tabla[i] = [xi,yi,K1,K2]
    if vertabla==True:
        np.set_printoptions(precision)
        titulo = ' [xi,     yi,     K1,    K2 ]'
        print(' EDO con Runge-Kutta 2do Orden primera derivada')
        print(titulo)
        print(tabla)
    return(tabla)

# PROGRAMA -----------------
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1
# INGRESO.
# 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

# PROCEDIMIENTO
tabla = rungekutta2(d1y,x0,y0,h,muestras,
                     vertabla=True)

# SALIDA
##print(' EDO con Runge-Kutta 2do Orden primera derivada')
##print(' [xi,     yi,     d1yi,',', K1,   K2 ]')
##print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
xi = tabla[:,0]
yi = tabla[:,1]
plt.plot(xi,yi)
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución Runge-Kutta 2do Orden primera derivada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
#plt.show() #comentar para la siguiente gráfica


# GRAFICA CON ANIMACION --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation

unmetodo = 'EDO: Runge-Kutta 2do Orden primera derivada'
narchivo = 'EdoRK2df' # nombre archivo GIF
muestras = 51 

# Puntos para la gráfica
a = xi[0]
b = xi[1]
K1 = tabla[:,2]
K2 = tabla[:,3]
n = len(xi)

# Parametros de trama/foto
retardo = 1000   # milisegundos entre tramas
tramas = len(xi)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
ymax = np.max([yi[0],yi[2]])
ymin = np.min([yi[0],yi[2]])
deltax = np.abs(xi[2]-xi[0])
deltay = np.abs(yi[2]-yi[0])
graf_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
# Lineas y puntos base
linea_fx, = graf_ani.plot(xi, yi,color='blue',
                          linestyle='dashed')
puntof, = graf_ani.plot(xi[0], yi[0],'o',
                        color='green',
                        label='xi,yi')
puntoa, = graf_ani.plot(xi[0], yi[0],'o',
                        color='Blue')
puntob, = graf_ani.plot(xi[1], yi[1],'o',
                        color='orange')

linea_h, = graf_ani.plot(xi, xi, color='orange',
                         label='h',
                         linestyle='dashed')
linea_K1, = graf_ani.plot(xi-0.02*deltax, xi-0.02*deltax,
                          color='green',label="K1",
                          linestyle='dashed')
linea_K2, = graf_ani.plot(xi+0.02*deltax, yi,
                          color='magenta',
                          label="K2",
                          linestyle='dashed')
linea_K12, = graf_ani.plot(xi, yi,
                          color='magenta')

# Cuadros de texto en gráfico
txt_i  = graf_ani.text(xi[0], yi[0]+0.05*deltay,'[x[i],y[i]]',
                       horizontalalignment='center')
txt_i1 = graf_ani.text(xi[1], xi[1]+0.05*deltay,'[x[i+1],y[i+1]]',
                       horizontalalignment='center')

# Configura gráfica
graf_ani.axhline(0, color='black')  # Linea horizontal en cero
graf_ani.set_title(unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Nueva Trama
def unatrama(i,xi,yi,term1,term2):   
    if i>1:
        ymax = np.max([yi[0:i+2]])
        ymin = np.min([yi[0:i+2]])
        deltax = np.abs(xi[i+1]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    else:
        ymax = np.max([yi[0],yi[2]])
        ymin = np.min([yi[0],yi[2]])
        deltax = np.abs(xi[2]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    # actualiza cada punto
    puntoa.set_xdata([xi[i]]) 
    puntoa.set_ydata([yi[i]])
    puntob.set_xdata([xi[i+1]]) 
    puntob.set_ydata([yi[i+1]])
    puntof.set_xdata(xi[0:i]) 
    puntof.set_ydata(yi[0:i])
    # actualiza cada linea
    linea_fx.set_xdata(xi[0:i+1])
    linea_fx.set_ydata(yi[0:i+1])
    linea_h.set_xdata([xi[i],xi[i+1]])
    linea_h.set_ydata([yi[i],yi[i]])
    linea_K1.set_xdata([xi[i+1]-0.02*deltax,xi[i+1]-0.02*deltax])
    linea_K1.set_ydata([yi[i],yi[i]+K1[i+1]])
    linea_K2.set_xdata([xi[i+1]+0.02*deltax,xi[i+1]+0.02*deltax])
    linea_K2.set_ydata([yi[i],yi[i]+K2[i+1]])
    linea_K12.set_xdata([xi[i+1]-0.02*deltax,xi[i+1]+0.02*deltax])
    linea_K12.set_ydata([yi[i]+K1[i+1],yi[i]+K2[i+1]])
    # actualiza texto
    txt_i.set_position([xi[i], yi[i]+0.05*deltay])
    txt_i1.set_position([xi[i+1], yi[i+1]+0.05*deltay])

    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_K1,linea_K2,linea_K12,
            txt_i,txt_i1,)
# Limpia Trama anterior
def limpiatrama(): 
    puntoa.set_ydata(np.ma.array(xi, mask=True))
    puntob.set_ydata(np.ma.array(xi, mask=True))
    puntof.set_ydata(np.ma.array(xi, mask=True))
    linea_h.set_ydata(np.ma.array(xi, mask=True))
    linea_K1.set_ydata(np.ma.array(xi, mask=True))
    linea_K2.set_ydata(np.ma.array(xi, mask=True))
    linea_K12.set_ydata(np.ma.array(xi, mask=True))
    
    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_K1,linea_K2,linea_K12,
            txt_i,txt_i1,)

# Trama contador
i = np.arange(0,tramas-1,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (xi,yi,K1,K2),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

    animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]
    ..

    Sistemas EDO. modelo depredador-presa con Runge-Kutta 2do Orden
    .

    Edo Presa Predador GIF animado

    Instrucciones en Python

    # Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision = 6):
    ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
        x0,y0 son valores iniciales, h es tamano de paso,
        muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+4),dtype=float)

    # incluye el punto [x0,y0,z0]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    i=0
    if vertabla==True:
        print('Runge-Kutta Segunda derivada')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        np.set_printoptions(precision)
        print(i,tabla[i,0:3])
        print('  ',tabla[i,3:])
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
        if vertabla==True:
            txt = ' '
            if i>=10:
                txt='  '
            print(str(i)+'',tabla[i,0:3])
            print(txt,tabla[i,3:])
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35

# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y

# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1

# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras,vertabla=True)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
print('Sistemas EDO: Modelo presa-predador')
##np.set_printoptions(precision=6)
##print(' [ ti, xi, yi]')
##print(tabla[:,0:4])

# GRAFICA tiempos vs población
import matplotlib.pyplot as plt

fig_t, (graf1,graf2) = plt.subplots(2)
fig_t.suptitle('Modelo predador-presa')
graf1.plot(ti,xi, color='blue',label='xi presa')

#graf1.set_xlabel('t tiempo')
graf1.set_ylabel('población x')
graf1.legend()
graf1.grid()

graf2.plot(ti,yi, color='orange',label='yi predador')
graf2.set_xlabel('t tiempo')
graf2.set_ylabel('población y')
graf2.legend()
graf2.grid()

# gráfica xi vs yi
fig_xy, graf3 = plt.subplots()
graf3.plot(xi,yi)

graf3.set_title('Modelo presa-predador [xi,yi]')
graf3.set_xlabel('x presa')
graf3.set_ylabel('y predador')
graf3.grid()
#plt.show()


# GRAFICA CON ANIMACION --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation
xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

unmetodo = 'Sistemas EDO Presa-Predador con Runge-Kutta'
narchivo = 'EdoPresaPredador' # nombre archivo GIF
muestras = 51 

# Puntos para la gráfica
a = xi[0]
b = xi[1]
n = len(ti)

# Parametros de trama/foto
retardo = 1000   # milisegundos entre tramas
tramas = len(xi)

# GRAFICA animada en fig_ani
fig_ani, (graf1_ani,graf2_ani) = plt.subplots(2)
ymax = np.max([yi[0],yi[2]])
ymin = np.min([yi[0],yi[2]])
deltax = np.abs(xi[2]-xi[0])
deltay = np.abs(yi[2]-yi[0])
graf1_ani.set_xlim([xi[0],xi[2]+0.05*deltax])
graf1_ani.set_ylim([ymin-0.05*deltay,ymax+0.05*deltay])

zmax = np.max([zi[0],zi[2]])
zmin = np.min([zi[0],zi[2]])
deltax = np.abs(xi[2]-xi[0])
deltaz = np.abs(zi[2]-zi[0])
graf2_ani.set_xlim([xi[0],xi[2]+0.05*deltax])
graf2_ani.set_ylim([zmin-0.05*deltaz,zmax+0.05*deltaz])
# Lineas y puntos base
linea_fx, = graf1_ani.plot(xi, yi,color='blue',
                          linestyle='dashed')
puntof, = graf1_ani.plot(xi[0], yi[0],'o',
                        color='blue',
                        label='xi,yi')
puntoa, = graf1_ani.plot(xi[0], yi[0],'o',
                        color='green')
puntob, = graf1_ani.plot(xi[1], yi[1],'o',
                        color='dodgerblue')
linea_h, = graf1_ani.plot(xi, xi, color='green',
                         label='h',
                         linestyle='dashed')
linea_term1, = graf1_ani.plot(xi, xi,
                             color='dodgerblue',label="(K1y+K2y)/2",
                             linestyle='dashed')
# Cuadros de texto en gráfico
#txt_i  = graf1_ani.text(xi[0], yi[0]+0.03*deltay,'[x[i],y[i]]',
#                       horizontalalignment='center')
#txt_i1 = graf1_ani.text(xi[1], xi[1]+0.03*deltay,'[x[i+1],y[i+1]]',
#                       horizontalalignment='center')

linea_gx, = graf2_ani.plot(xi, zi,color='orange',
                          linestyle='dashed')
puntog, = graf2_ani.plot(xi[0], zi[0],'o',
                        color='orange',
                        label='xi,zi')
puntog_a, = graf2_ani.plot(xi[0], zi[0],'o',
                        color='green')
puntog_b, = graf2_ani.plot(xi[1], zi[1],'o',
                        color='red')
lineag_h, = graf2_ani.plot(xi, xi, color='green',
                         label='h',
                         linestyle='dashed')
lineag_term1, = graf2_ani.plot(xi, xi,
                             color='red',label="(K1z+K2z)/2",
                             linestyle='dashed')

# Configura gráfica
graf1_ani.axhline(0, color='black')  # Linea horizontal en cero
graf1_ani.set_title(unmetodo)
graf1_ani.set_xlabel('x')
graf1_ani.set_ylabel('y(x)')
graf1_ani.legend()
graf1_ani.grid()

graf2_ani.axhline(0, color='black')  # Linea horizontal en cero
#graf2_ani.set_title(unmetodo)
graf2_ani.set_xlabel('x')
graf2_ani.set_ylabel('z(x)')
graf2_ani.legend()
graf2_ani.grid()

# Nueva Trama
def unatrama(i,xi,yi,zi):
    
    if i>1:
        ymax = np.max([yi[0:i+2]])
        ymin = np.min([yi[0:i+2]])
        deltax = np.abs(xi[i+1]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf1_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf1_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])

        zmax = np.max([zi[0:i+2]])
        zmin = np.min([zi[0:i+2]])
        deltaz = np.abs(zmax-zmin)
        graf2_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf2_ani.set_ylim([zmin-0.05*deltaz,zmax+0.1*deltaz])
    else:
        ymax = np.max([yi[0:2]])
        ymin = np.min([yi[0:2]])
        deltax = np.abs(xi[2]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf1_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf1_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])

        zmax = np.max([zi[0:2]])
        zmin = np.min([zi[0:2]])
        deltaz = np.abs(zmax-zmin)
        graf2_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf2_ani.set_ylim([zmin-0.05*deltaz,zmax+0.1*deltaz])
    # actualiza cada punto
    puntoa.set_xdata([xi[i]]) 
    puntoa.set_ydata([yi[i]])
    puntob.set_xdata([xi[i+1]]) 
    puntob.set_ydata([yi[i+1]])
    puntof.set_xdata(xi[0:i]) 
    puntof.set_ydata(yi[0:i])
    # actualiza cada linea
    linea_fx.set_xdata(xi[0:i+1])
    linea_fx.set_ydata(yi[0:i+1])
    linea_h.set_xdata([xi[i],xi[i+1]])
    linea_h.set_ydata([yi[i],yi[i]])
    linea_term1.set_xdata([xi[i+1],xi[i+1]])
    linea_term1.set_ydata([yi[i],yi[i+1]])
    # actualiza texto
    #txt_i.set_position([xi[i], yi[i]+0.03*deltay])
    #txt_i1.set_position([xi[i+1], yi[i+1]+0.03*deltay])

    # actualiza cada punto
    puntog_a.set_xdata([xi[i]]) 
    puntog_a.set_ydata([zi[i]])
    puntog_b.set_xdata([xi[i+1]]) 
    puntog_b.set_ydata([zi[i+1]])
    puntog.set_xdata(xi[0:i]) 
    puntog.set_ydata(zi[0:i])
    # actualiza cada linea
    linea_gx.set_xdata(xi[0:i+1])
    linea_gx.set_ydata(zi[0:i+1])
    lineag_h.set_xdata([xi[i],xi[i+1]])
    lineag_h.set_ydata([zi[i],zi[i]])
    lineag_term1.set_xdata([xi[i+1],xi[i+1]])
    lineag_term1.set_ydata([zi[i],zi[i+1]])

    return (puntoa,puntob,puntof,
            linea_fx,linea_h,
            linea_term1,)
            #txt_i,txt_i1,)
# Limpia Trama anterior
def limpiatrama(): 
    puntoa.set_ydata(np.ma.array(xi, mask=True))
    puntob.set_ydata(np.ma.array(xi, mask=True))
    puntof.set_ydata(np.ma.array(xi, mask=True))
    linea_h.set_ydata(np.ma.array(xi, mask=True))
    linea_term1.set_ydata(np.ma.array(xi, mask=True))

    puntog_a.set_ydata(np.ma.array(xi, mask=True))
    puntog_b.set_ydata(np.ma.array(xi, mask=True))
    puntog.set_ydata(np.ma.array(xi, mask=True))
    lineag_h.set_ydata(np.ma.array(xi, mask=True))
    lineag_term1.set_ydata(np.ma.array(xi, mask=True))
    return (puntoa,puntob,puntof,
            linea_fx,linea_h,
            linea_term1,
            puntog_a,puntog_b,puntog,
            linea_gx,lineag_h,
            lineag_term1,)
            #txt_i,txt_i1,)

# Trama contador
i = np.arange(0,tramas-1,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (xi,yi,zi),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

    animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

  • 6.4 Sistemas EDO. modelo depredador-presa con Runge-Kutta y Python

    Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ]
    ..

    1. Ejercicio

    Referencia: Chapra 28.2 p831 pdf855, Rodríguez 9.2.1 p263
    https://es.wikipedia.org/wiki/Ecuaciones_Lotka%E2%80%93Volterra
    https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations

    Modelos depredador-presa y caos. Ecuaciones Lotka-Volterra. En el sistema de ecuaciones:

    \frac{dx}{dt} = ax - bxy \frac{dy}{dt} = -cy + dxy

    predador Presa 01variables
    x = número de presas
    y = número de depredadores
    t = tiempo de observación
    coeficientes
    a = razón de crecimiento de la presa, (0.5)
    c = razón de muerte del depredador (0.35)
    b = efecto de la interacción depredador-presa sobre la muerte de la presa (0.7)
    d = efecto de la interacción depredador-presa sobre el crecimiento del depredador, (0.35)

    Considere como puntos iniciales en la observación de las especies:
    t = 0, x = 2, y = 1, h = 0.5

    Los términos que multiplican xy hacen que las ecuaciones sean no lineales.
    Observe que la variable tiempo no se encuentra en las expresiones f y g, h se aplica a tiempo.

    Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ]
    ..

    2. Desarrollo analíticoEdo Presa Predador GIF animado

    Para resolver el sistema, se plantean las ecuaciones de forma simplificada para el algoritmo:

    f(t,x,y) = 0.5 x - 0.7 xy g(t,x,y) = -0.35y + 0.35xy

    Las expresiones se adaptan al método de Runge-Kutta para primeras derivadas por cada variable de población. Se deben usar de forma simultánea para cada tiempo t.

    K1x = h f(t,x,y) = 0.5 \Big( 0.5 x - 0.7 xy \Big) K1y = h g(t,x,y) = 0.5 \Big(-0.35y + 0.35xy \Big)

    ..

    K2x = h f(t+h,x+K1x,y+K1y) = 0.5 \Big( 0.5 (x+K1x) - 0.7 (x+K1x)(y+K1y) \Big) K2y = h g(t+h,x+K1x,y+K1y) = 0.5 \Big(-0.35(y+K1y) + 0.35(x+K1x)(y+K1y) \Big)

    ..

    x[i+1] = x[i] + \frac{K1x+K2x}{2} y[i+1] = y[i] + \frac{K1y+K2y}{2} t[i+1] = t[i] + h

    con lo que se puede aplicar al ejercicio en cada iteración. dadas las condiciones iniciales.

    itera = 0

    t = 0, x = 2, y = 1, h = 0.5

    K1x = 0.5 f(0,2,1) = 0.5 \Big( 0.5 (2) - 0.7 (2)(1) \Big) = -0.2 K1y = 0.5 g(0,2,1) = 0.5 \Big(-0.35(1) + 0.35(2)(1) \Big) =0.175

    ..

    K2x = 0.5 f(0+0.5, 2+(-0.2), 1+0.175) = 0.5 \Big( 0.5 (2+(-0.2)) - 0.7 (2+(-0.2))(1+0.175) \Big) = -0.29025 K2y = 0.5 g(0+0.5, 2+(-0.2), 1+0.175) = 0.5 \Big(-0.35(1+0.175) + 0.35(2+(-0.2))(1+0.175) \Big) = 0.1645

    ..

    x[1] = x[0] + \frac{K1x+K2x}{2} = 2 + \frac{-0.2+(-0.29025)}{2} = 1.7548 y[1] = y[0] + \frac{K1y+K2y}{2} = 1 + \frac{0.175+0.1645}{2}= 1.1697 t[1] = t[0] + h = 0 +0.5 = 0.5

    itera = 1

    t = 0.5, x = 1.7548, y = 1.1697, h = 0.5

    K1x = 0.5 \Big( 0.5 (0,1.7548) - 0.7 (0,1.7548)(1.1697) \Big) = -0.2797 K1y = 0.5 \Big(-0.35(1.1697) + 0.35(1.7548)(1.1697) \Big) =0.1545

    ..

    K2x = 0.5 \Big( 0.5 (1.7548+(-0.2797)) - 0.7 (1.7548+(-0.2797))(1.1697+0.1545) \Big) =-0.3149 K2y = 0.5 \Big(-0.35(1.1697+0.1545) + 0.35(1.7548+(-0.2797))(1.1697+0.1545) \Big) = 0.1645

    ..

    x[2] = 1.7548 + \frac{-0.2797+(-0.3149)}{2} = 1.4575 y[2] = 1.1697 + \frac{0.1545+0.1645}{2} = 1.3020 t[2] = t[0] + h = 0.5 +0.5 = 1

    itera=2

    t = 1, x = 1.4575, y = 1.3020, h = 0.5

    continuar como tarea ...

    Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ]

    ..

    3. Algoritmo en Python

    Planteamiento que se ingresan al algoritmo con el algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en

    EDO con Runge-Kutta d2y/dx2

    Al ejecutar el algoritmo se obtienen los siguientes resultados:

    Runge-Kutta Segunda derivada
i  [ xi,  yi,  zi ]
   [ K1y,  K1z,  K2y,  K2z ]
0 [0. 2. 1.]
   [0. 0. 0. 0.]
1 [0.5      1.754875 1.16975 ]
  [-0.2      0.175   -0.29025  0.1645 ]
2 [1.       1.457533 1.302069]
  [-0.279749  0.154528 -0.314935  0.11011 ]
3 [1.5      1.167405 1.373599]
  [-0.29985   0.104254 -0.280406  0.038807]
4 [2.       0.922773 1.381103]
  [-0.26939   0.040241 -0.219874 -0.025233]
5 [2.5      0.734853 1.33689 ]
  [-0.215362 -0.018665 -0.160478 -0.069761]
6 [3.       0.598406 1.258434]
  [-0.160133 -0.062033 -0.11276  -0.09488 ]
...

    Los resultados de la tabla se muestran parcialmente, pues se usaron mas de 100 iteraciones.

    Los resultados se pueden observar de diferentes formas:

    a) Cada variable xi, yi versus ti, es decir cantidad de animales de cada especie durante el tiempo de observación

    predador Presa 02

    b) Independiente de la unidad de tiempo, xi vs yi, muestra la relación entre la cantidad de presas y predadores. Relación que es cíclica y da la forma a la gráfica.

    predador Presa 03

    Las instrucciones del algoritmo en Python usadas en el problema son:

    # Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision=6):
    ''' solucion a EDO d2y/dx2 con Runge-Kutta 2do Orden,
    f(x,y,z) = z #= y'
    g(x,y,z) = expresion d2y/dx2 con z=y'
    tambien es solucion a sistemas edo f() y g()
    x0,y0,z0 son valores iniciales, h es tamano de paso,
    muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+4),dtype=float)
    # incluye el punto [x0,y0,z0,K1y,K1z,K2y,K2z]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    
    xi = x0 # valores iniciales
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
        
    if vertabla==True:
        np.set_printoptions(precision)
        print('EDO f,g con Runge-Kutta 2 Orden')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        for i in range(0,tamano,1):  
            txt = ' '
            if i>=10:
                txt = '  '
            print(str(i),tabla[i,0:3])
            print(txt,tabla[i,3:])
    
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35

# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y

# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1

# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras,vertabla=True)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
print('Sistemas EDO: Modelo presa-predador')
##np.set_printoptions(precision=6)
##print(' [ ti, xi, yi]')
##print(tabla[:,0:4])

    Los resultados numéricos se usan para generar las gráficas presentadas, añadiendo las instrucciones:

    # GRAFICA tiempos vs población
import matplotlib.pyplot as plt

fig_t, (graf1,graf2) = plt.subplots(2)
fig_t.suptitle('Modelo predador-presa')
graf1.plot(ti,xi, color='blue',label='xi presa')

#graf1.set_xlabel('t tiempo')
graf1.set_ylabel('población x')
graf1.legend()
graf1.grid()

graf2.plot(ti,yi, color='orange',label='yi predador')
graf2.set_xlabel('t tiempo')
graf2.set_ylabel('población y')
graf2.legend()
graf2.grid()

# gráfica xi vs yi
fig_xy, graf3 = plt.subplots()
graf3.plot(xi,yi)

graf3.set_title('Modelo presa-predador [xi,yi]')
graf3.set_xlabel('x presa')
graf3.set_ylabel('y predador')
graf3.grid()
plt.show()

    Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ]

  • 6.3 EDO d2y/dx2, Runge-Kutta con Python

    [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
    Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
    ..

    1. EDO \frac{\delta^2 y}{\delta x^2} con Runge-Kutta

    Para una ecuación diferencial de segunda derivada (segundo orden) con condiciones de inicio en x0, y0, y'0

    \frac{\delta ^2 y}{\delta x^2} = \frac{\delta y}{\delta x} + etc

    Forma estandarizada de la ecuación:

    y'' = y' + etc

    Se puede sustituir la variable y' por z, lo que se convierte a dos expresiones que forman un sistema de ecuaciones:

    \begin{cases} z= y' = f_x(x,y,z) \\ z' = (y')' = z + etc = g_x(x,y,z) \end{cases}

    y se pueden reutilizar los métodos para primeras derivadas, por ejemplo Runge-Kutta de 2do y 4to orden para las variables x,y,z de forma simultanea.

    Runge-Kutta 2do Orden tiene error de truncamiento O(h3)
    Runge-Kutta 4do Orden tiene error de truncamiento O(h5)

    Al aplicar Runge-Kutta para las variables dependientes, se tiene:

    y'' = y' + etc \begin{cases} f_x(x,y,z) = z \\ g_x(x,y,z) = z + etc \end{cases} K_{1y} = h f(x_i, y_i, z_i) K_{1z} = hg(x_i, y_i, z_i) K_{2y} = h f(x_i +h, y_i + K_{1y} , z_i + K_{1z}) K_{2z} = h g(x_i +h, y_i + K_{1y}, z_i + K_{1z}) y_{i+1}=y_i+\frac{K_{1y}+K_{2y}}{2} z_{i+1}=z_i+\frac{K_{1z}+K_{2z}}{2} x_{i+1} = x_i +h

    [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
    Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]

    ..

    2. Ejercicio

    Referencia: Chapra Ejercicio 25.23 p265, 2Eva_IT2018_T1 Paracaidista wingsuit http://www.elperiodicodearagon.com/noticias/sociedad/alarma-francia-cinco-muertes-verano-moda-hombres-pajaro-wingsuit_877164.html

    Si suponemos que el arrastre es proporcional al cuadrado de la velocidad, se puede modelar la altura (y) de un objeto que cae, como un paracaidista, por medio de la ecuación diferencial ordinaria:

    \frac{d^2y}{dt^2} = g - \frac{Cd}{m} \Big( \frac{dy}{dt} \Big) ^2

    Que es una EDO de 2do orden o como 2da derivada.

    Resuelva para la altura que recorre un objeto de 90 Kg con coeficiente de arrastre  Cd =0.225 kg/m.

    Si la velocidad inicial es 0 y la altura inicial es 1 Km, determine la velocidad y posición en cada tiempo, usando un tamaño de paso de 2s.

    [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
    Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]

    ..

    3. Desarrollo analítico

    La solución propone resolver de forma simultanea para t,y,v con Runge Kutta para segunda derivada donde:diagrama cuerpo libre

    f(t,y,v) = -v g(t,y,v) = g - \frac{Cd}{m} v^2

    Al sustituir los valores de las constantes en la ecuación como gravedad, masa e índice de arrastre se tiene:

    f(t,y,v) = -v g(t,y,v) = 9.8 - \frac{0.225}{90} v^2

    con las condiciones iniciales del ejercicio  t0 = 0 , y0 = 1000, v0 = 0
    la velocidad se inicia con cero, si el paracaidista se deja caer desde el borde el risco, como en el video adjunto al enunciado.

    Para las iteraciones, recuerde que
    t se corrige con t+h (en el algoritmo era la posición para x)
    y se corrige con y+K1y
    v se corrige con v+K1v (en el algoritmo era la posición para z)

    itera = 0

    K1y = h f(t,y,v) = 2(-(0)) = 0 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (0)^2) = 19.6

    ..
    K2y = h f(t+h, y+K1y, v + K1v) = 2(-(0 + 19.6)) = -39.2

    K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (0+19.6)^2) =17.6792

    ..
    y_1 = y_0 + \frac{K1y+K2y}{2} = 1000 + \frac{0-39.2}{2}= 980.4

    v_1 = v_0 + \frac{K1v+K2v}{2} = 0 + \frac{19.6-17.67}{2} = 18.63 t_1 =t_0 + h = 0+2 = 2
    ti yi vi K1y K1v K2y K2v
    0 1000 0 0 0 0 0
    2 980.4 18.63 0 19.6 -39.2 17.6792

    itera = 1

    K1y = 2(-(18.63)) = -37.2792 K1v = 2(9.8 - \frac{0.225}{90} (18.63)^2) = 17.8628

    ..
    K2y =2(-(18.6396+17.8628)) =-73.00

    K2v = 2(9.8 - \frac{0.225}{90} (18.6396+17.8628)^2) =12.9378

    ..
    y_2 =980.4 + \frac{ -37.2792+(-73.00)}{2}= 925.25

    v_2 = 18.63 + \frac{17.8628+12.9378}{2} = 34.0399 t_2 =t_1 + h = 2+2 = 4
    ti yi vi K1y K1v K2y K2v
    0 1000 0 0 0 0 0
    2 980.4 18.63 0 19.6 -39.2 17.6792
    4 925.25 34.0399 -37.2792 17.8628 -73.00 12.9378

    itera = 2

    K1y = h f(t,y,v) = 2(-(34.0399)) = -68.0798 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (34.0399)^2) = 13.8064

    ..
    K2y = h f(t+h, y+K1y, v + K1v) = 2(-(34.0399+13.8064)) =-95.6927

    K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (34.0399+13.8064)^2) =8.1536

    ..
    y_2 = 925.25 + \frac{ -68.0798+(-95.6927)}{2}= 843.3716

    v_2 = 34.0399 + \frac{13.8064+8.1536}{2} = 45.0199 t_2 =t_1 + h = 4+2 = 6
    ti yi vi K1y K1v K2y K2v
    0 1000 0 0 0 0 0
    2 980.4 18.63 0 19.6 -39.2 17.6792
    4 925.25 34.0399 -37.2792 17.8628 -73.00 12.9378
    6 843.3716 45.0199 -68.0798 13.8064 -95.6927 8.1536

    Usando el algoritmo se obtiene el siguiente resultado:

    EDO f,g con Runge-Kutta 2 Orden
i  [ xi,  yi,  zi ]
   [ K1y,  K1z,  K2y,  K2z ]
0 [   0. 1000.    0.]
  [0. 0. 0. 0.]
1 [  2.     980.4     18.6396]
  [  0.      19.6    -39.2     17.6792]
2 [  4.     925.258   34.0399]
  [-37.2792  17.8628 -73.0049  12.9379]
3 [  6.     843.3717  45.02  ]
  [-68.0799  13.8064 -95.6927   8.1536]
4 [  8.     743.8657  52.1312]
  [ -90.0399    9.466  -108.972     4.7564]
5 [ 10.     633.5917  56.4855]
  [-104.2623    6.0117 -116.2857    2.697 ]
6 [ 12.     516.9737  59.0692]

    con la siguiente gráfica

    paracaidista wingsuit grafica 02

    [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
    Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]

    ..

    4. Algoritmo en Python para  \frac{\delta^2 y}{\delta x^2} Runge-Kutta 2do Orden

    Se presenta las instrucciones en Python con una función.

    # EDO d2y/dx2. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0,z0 entrega tabla[xi,yi,zi,K1y,K1z,K2y,K2z]
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision=6):
    ''' solucion a EDO d2y/dx2 con Runge-Kutta 2do Orden,
    f(x,y,z) = z #= y'
    g(x,y,z) = expresion d2y/dx2 con z=y'
    tambien es solucion a sistemas edo f() y g()
    x0,y0,z0 son valores iniciales, h es tamano de paso,
    muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+4),dtype=float)
    # incluye el punto [x0,y0,z0,K1y,K1z,K2y,K2z]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    
    xi = x0 # valores iniciales
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
        
    if vertabla==True:
        np.set_printoptions(precision)
        print('EDO f,g con Runge-Kutta 2 Orden')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        for i in range(0,tamano,1):  
            txt = ' '
            if i>=10:
                txt = '  '
            print(str(i),tabla[i,0:3])
            print(txt,tabla[i,3:])
    
    return(tabla)

# PROGRAMA PRUEBA -------------------
# 2Eva_IT2018_T1 Paracaidista wingsuit

# INGRESO
f = lambda t,y,v: -v # el signo, revisar diagrama cuerpo libre
g = lambda t,y,v: 9.8 - (0.225/90)*(v**2)
t0 = 0
y0 = 1000
v0 = 0
h  = 2
muestras = 15+1

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras,
                       vertabla=True, precision=4)
ti = tabla[:,0]
yi = tabla[:,1]
vi = tabla[:,2]

# SALIDA
# print('tabla de resultados')
# print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
plt.subplot(211)
plt.plot(ti,vi,label='velocidad v(t)', color='green')
plt.plot(ti,vi,'o')
plt.ylabel('velocidad (m/s)')
plt.title('paracaidista Wingsuit con Runge-Kutta')
plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(ti,yi,label='Altura y(t)',)
plt.plot(ti,yi,'o',)
plt.axhline(0, color='red')
plt.xlabel('tiempo (s)')
plt.ylabel('Altura (m)')
plt.legend()
plt.grid()

plt.show()

    [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
    Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
    ..

    5. Algoritmo en Python para  \frac{\delta^2 y}{\delta x^2} Runge-Kutta 4to Orden

    Se adjunta la función para 4to orden que se puede probar con el mismo ejercicio anterior.

    # EDO d2y/dx2. Método de RungeKutta 4to Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0,z0 entrega
# tabla[xi,yi,zi,K1y,K1z,K2y,K2z,K3y,K3z,K4y,K4z]
import numpy as np

def rungekutta4_fg(fx,gx,x0,y0,z0,h,muestras,
                   vertabla=False, precision=6):
    ''' solucion a EDO d2y/dx2 con Runge-Kutta 4to Orden,
    f(x,y,z) = z #= y'
    g(x,y,z) = expresion d2y/dx2 con z=y'
    tambien es solucion a sistemas edo f() y g()
    x0,y0,z0 son valores iniciales, h es tamano de paso,
    muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+8),dtype=float)
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,z0,0,0,0,0,0,0,0,0]

    xi = x0 # valores iniciales
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
        K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
        
        K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
        K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)

        K4y = h * fx(xi+h, yi + K3y, zi + K3z)
        K4z = h * gx(xi+h, yi + K3y, zi + K3z)

        yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
        zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
        xi = xi + h
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z,K3y,K3z,K4y,K4z]
    
    if vertabla==True:
        np.set_printoptions(precision)
        print('EDO f,g con Runge-Kutta 4 Orden')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        print('   [ K3y,  K3z,  K4y,  K4z ]')
        for i in range(0,tamano,1):  
            txt = ' '
            if i>=10:
                txt = '  '
            print(str(i),tabla[i,0:3])
            print(txt,tabla[i,3:7])
            print(txt,tabla[i,7:])

    return(tabla)

    [ Runge Kutta  \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
    Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]

  • 6.2.1 EDO dy/dx, Runge-Kutta 4to Orden con Python

    [ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

    ..

    EDO  \frac{\delta y}{\delta x} Runge-Kutta 4to Orden

    Referencia: Chapra 25.3.3 p746, Rodríguez 9.1.8 p358

    Para una ecuación diferencial de primera derivada (primer orden) con una condición de inicio:
    Runge Kutta 4to Orden

    \frac{\delta y}{\delta x} + etc =0 y'(x) = f(x_i,y_i) y(x_0) = y_0

    La fórmula de Runge-Kutta de 4to orden realiza una corrección con 4 valores de K:

    y_{i+1} = y_i + \frac{K_1 + 2K_2 + 2K_3 + K_4}{6}

    debe ser equivalente a la serie de Taylor de 5 términos:

    y_{i+1} = y_i + h f(x_i,y_i) + + \frac{h^2}{2!} f'(x_i,y_i) + \frac{h^3}{3!} f''(x_i,y_i) + +\frac{h^4}{4!} f'''(x_i,y_i) + O(h^5) x_{i+1} = x_i + h

    Runge-Kutta 4do Orden tiene error de truncamiento O(h5)

    Ejercicio

    Para el desarrollo analítico se tienen las siguientes expresiones para el ejercicio usado en Runge-Kutta de orden 2, que ahora será con orden 4:

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

    Se usa las expresiones de Runge-Kutta en orden, K1 corresponde a una corrección de EDO con Taylor de dos términos (método de Euler). K2 considera el cálculo a medio tamaño de paso más adelante.

    iteración:

    K_1 = h f(x_i,y_i) = 0.1 (y_i -x_i^2 +x_i +1) K_2 = h f\Big(x_i+\frac{h}{2}, y_i + \frac{K_1}{2} \Big) K_2 = 0.1 \Big(\big(y_i+\frac{K_1}{2}\big) -\big(x_i+\frac{h}{2}\big)^2 +\big(x_i+\frac{h}{2}\big) +1 \Big) K_3 = h f\Big(x_i+\frac{h}{2}, y_i + \frac{K_2}{2} \Big) K_3 = 0.1 \Big(\big(y_i+\frac{K_2}{2}\big) -\big(x_i+\frac{h}{2}\big)^2 +\big(x_i+\frac{h}{2}\big) +1 \Big) K_4 = h f(x_i+h, y_i + K_3 ) K_4 = 0.1 \Big((y_i+K_3) -(x_i+h)^2 +(x_i+h) +1 \Big) y_{i+1} = y_i + \frac{K_1+2K_2+2K_3+K_4}{6} x_{i+1} = x_i + h

    Las iteraciones se dejan como tarea

    [ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

    ..

    Algoritmo en Python como Función

    def rungekutta4(d1y,x0,y0,h,muestras, vertabla=False, precision=6):
    ''' solucion a EDO con Runge-Kutta 4do Orden primera derivada,
        x0,y0 son valores iniciales, tamaño de paso h.
        muestras es la cantidad de puntos a calcular.
    '''
    # Runge Kutta de 4do orden
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+4),dtype=float)
    
    # incluye el punto [x0,y0,K1,K2,K3,K4]
    tabla[0] = [x0,y0,0,0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h/2, yi + K1/2)
        K3 = h * d1y(xi+h/2, yi + K2/2)
        K4 = h * d1y(xi+h, yi + K3)

        yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
        xi = xi + h
        
        tabla[i] = [xi,yi,K1,K2,K3,K4]
        
    if vertabla==True:
        np.set_printoptions(precision)
        titulo = ' [xi,     yi,     K1,    K2,     K3,     K4 ]'
        print(' EDO con Runge-Kutta 4do Orden primera derivada')
        print(titulo)
        print(tabla)
    return(tabla)

    Note que el método de Runge-Kutta de 4to orden es similar a la regla de Simpson 1/3. La ecuación representa un promedio ponderado para establecer la mejor pendiente.

    [ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

    ..

    Ejercicio

    2Eva_IT2018_T1 Paracaidista wingsuit

    Solución Propuesta: s2Eva_IT2018_T1 Paracaidista wingsuit

     

    La segunda parte corresponde a Runge-Kutta de 4to Orden

    [ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]