7.2.1 EDP Elípticas método iterativo

con el resultado desarrollado en EDP elípticas para:

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 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 numerico 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.

# Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 60
Tb = 60
Tc = 50
Td = 70
# dimensiones de la placa
x0 = 0
xn = 2
y0 = 0
yn = 1.5
# discretiza, supone dx=dy
dx = 0.25 
dy = 0.25 
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)
# Matriz u
u = np.zeros(shape=(n,m),dtype = float)
# valores en fronteras
u[0,:]   = Ta
u[n-1,:] = Tb
u[:,0]   = Tc
u[:,m-1] = Td

# valor inicial de iteración
promedio = (Ta+Tb+Tc+Td)/4
u[1:n-1,1:m-1] = promedio
# iterar puntos interiores
itera = 0
converge = 0
while not(itera>=maxitera or converge==1):
    itera = itera +1
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1])/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge = 1

# SALIDA
np.set_printoptions(precision=2)
print('converge = ', converge)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

Un ejemplo de resultados:

converge = 1
xi=
[ 0.    0.25  0.5   0.75  1.    1.25  1.5   1.75  2.  ]
yi=
[ 0.    0.25  0.5   0.75  1.    1.25  1.5 ]
matriz u
[[ 50.    60.    60.    60.    60.    60.    70.  ]
 [ 50.    55.6   58.23  60.    61.77  64.4   70.  ]
 [ 50.    54.15  57.34  60.    62.66  65.85  70.  ]
 [ 50.    53.67  56.97  60.    63.03  66.33  70.  ]
 [ 50.    53.55  56.87  60.    63.13  66.45  70.  ]
 [ 50.    53.67  56.97  60.    63.03  66.33  70.  ]
 [ 50.    54.15  57.34  60.    62.66  65.85  70.  ]
 [ 50.    55.6   58.23  60.    61.77  64.4   70.  ]
 [ 50.    60.    60.    60.    60.    60.    70.  ]]
>>>

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

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 graficada se obtiene como la transpuesta de u

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x

figura = plt.figure()
ax = Axes3D(figura)
ax.plot_surface(X, Y, U,
                rstride=1,
                cstride=1,
                cmap=cm.Reds)

plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

7.1.3 EDP Parabólicas método implícito

Referencia:  Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodriguez 10.2.4 p415

Se utiliza el mismo ejercicio presentado en método explícito.

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

En éste caso se usan diferencias finitas centradas y hacia atras; 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.

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.

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.


Algoritmos en Python

para la solución con el método implícito.

# 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
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 60
Tb = 40
T0 = 25
# longitud en x
a = 0
b = 1
# Constante K
K = 4
# Tamaño de paso
dx = 0.1
dt = 0.01
# iteraciones en tiempo
n = 100

# PROCEDIMIENTO
# Valores de x
xi = np.arange(a,b+dx,dx)
m  = len(xi)
ultimox = m-1

# Resultados en tabla de u
u = np.zeros(shape=(m,n), dtype=float)
# valores iniciales de u[:,j]
j = 0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0
u[ultimox,j] = Tb

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

# Calcula U para cada tiempo + dt
j = 1
while not(j>=n):
    u[0,j] = Ta
    u[m-1,j] = Tb

    # Matriz de ecuaciones
    tamano = m-2
    A = np.zeros(shape=(tamano,tamano), dtype = float)
    B = np.zeros(tamano, dtype = float)
    for f in range(0,tamano,1):
        if (f>0):
            A[f,f-1]=P
        A[f,f] = Q
        if (f<(tamano-1)):
            A[f,f+1]=R
        B[f] = -u[f+1,j-1]
    B[0] = B[0]-P*u[0,j]
    B[tamano-1] = B[tamano-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,tamano,1):
        u[f+1,j] = C[f]

    # siguiente iteración
    j = j + 1
        
# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

algunos valores:

Tabla de resultados
[[ 60.    60.    60.   ...,  60.    60.    60.  ]
 [ 25.    31.01  35.25 ...,  57.06  57.09  57.11]
 [ 25.    26.03  27.49 ...,  54.22  54.26  54.31]
 ..., 
 [ 25.    25.44  26.07 ...,  42.22  42.27  42.31]
 [ 25.    27.57  29.39 ...,  41.07  41.09  41.11]
 [ 40.    40.    40.   ...,  40.    40.    40.  ]]

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

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.m')

plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
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 pdf915, Burden 9Ed 12.2 p727, Rodriguez 10.2.2 pdf409 .

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:

http://blog.espol.edu.ec/analisisnumerico/recursos/resumen-python/tiempos-de-ejecucion/

7.1.1 EDP Parabólicas método explícito

Referencia:  Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodriguez 10.2 p406

Siguiendo el tema anterior en 9.1, se tiene que:

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

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

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.

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, Rodriguez 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 suponga:

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

 


Algoritmo en Python

Se presenta la propuesta en algoritmo para el método explícito.

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
# Referencia: Chapra 30.2 p.888 pdf.912
#       Rodriguez 10.2 p.406
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 60
Tb = 40
T0 = 25
# longitud en x
a = 0
b = 1
# Constante K
K = 4
# Tamaño de paso
dx = 0.1
dt = dx/10
# iteraciones en tiempo
n = 200

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

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,:]= Ta
u[1:ultimox,j] = T0
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
j = 0
while not(j>=ultimot): # igual con lazo for
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

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

Tabla de resultados
[[ 60.    60.    60.   ...,  60.    60.    60.  ]
 [ 25.    33.75  38.12 ...,  57.93  57.93  57.93]
 [ 25.    25.    27.19 ...,  55.86  55.86  55.87]
 ..., 
 [ 25.    25.    25.94 ...,  43.86  43.86  43.87]
 [ 25.    28.75  30.62 ...,  41.93  41.93  41.93]
 [ 40.    40.    40.   ...,  40.    40.    40.  ]]
>>> 

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.

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
    
plt.xlabel('x[i]')
plt.ylabel('t[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 segúndo índice es columna.

Tarea o Proyecto: Realizar la animación de los cambios de temperatura en el tiempo.

6.3 Sistemas EDO. modelo depredador-presa

Referencia: Chapra 28.2 p831 pdf855, Rodriguez 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

variables
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)

Los términos que multiplican xy hacen que las ecuaciones sean no lineales.

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

f =  a x - b x y
g = -c y + d x y

con valores iniciales: t = 0, x = 2, y = 1, h = 0.5

Observe que la variable tiempo no se encuentra en las expresiones f y g, h se aplica a tiempo.


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

Runge-Kutta d2y/dx2

Al ejecutar el algoritmo se obtienen los siguientes resultados:

 [ ti, xi, yi] 
[[  0.         2.         1.      ]
 [  0.5        1.754875   1.16975 ]
 [  1.         1.457533   1.302069]
 [  1.5        1.167405   1.373599]
 [  2.         0.922773   1.381103]
 [  2.5        0.734853   1.33689 ]
 [  3.         0.598406   1.258434]
 ...
 [ 49.         1.11309    1.389894]
 [ 49.5        0.876914   1.38503 ]
 [ 50.         0.698717   1.331107]
 [ 50.5        0.570884   1.246132]]
>>> 

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

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.


Algoritmo con Python

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,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,3),dtype=float)
    tabla[0] = [t0,x0,y0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi]
    tabla = np.array(tabla)
    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)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi]')
print(tabla)

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

plt.plot(ti,xi, label='xi presa')
plt.plot(ti,yi, label='yi predador')

plt.title('Modelo predador-presa')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

# gráfica xi vs yi
plt.plot(xi,yi)

plt.title('Modelo presa-predador [xi,yi]')
plt.xlabel('x presa')
plt.ylabel('y predador')
plt.grid()
plt.show()

Tarea: Añadir la animación de la gráfica usando la variable tiempo para mostrar un punto que marque el par ordenado[x,y] al variar t.

6.2.2 EDO Runge-Kutta d2y/dx2

Para sistemas de ecuaciones diferenciales ordinarias 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 convierte a:

z= y' = f_x(x,y,z) z' = (y')' = z + etc = g_x(x,y,z)

con las condiciones de inicio en x0, y0, z0

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

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

Runge-Kutta 2do Orden d2y/dx2 en Python:

import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)

    # incluye el punto [x0,y0,z0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    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
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

Runge-Kutta 4do Orden d2y/dx2 en Python:

def rungekutta4_fg(fx,gx,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)

    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    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
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

Una aplicación del algoritmo en Señales y Sistemas:

LTI CT – Respuesta entrada cero – Desarrollo analítico, TELG1001-Señales y Sistemas

6.2.1 EDO Runge-Kutta 4to Orden dy/dx

Referencia: Chapra 25.3.3 p746 pdf 770, Rodriguez 9.1.8 p358

Para una ecuación diferencial de primer orden con una condición de inicio, la fórmula de Runge-Kutta de 4to orden se obtiene de la expresión con cinco términos:

y_{i+1} = y_i + aK_1 + bK_2 + cK_3 + dK_4

siendo:

y'(x) = f(x_i,y_i) y(x_0) = y_0

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

que usando aproximaciones de derivadas, se obtienen:

# Runge Kutta de 4do orden
def rungekutta4(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    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
        
        estimado[i] = [xi,yi]
    return(estimado)

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.

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

6.2 EDO Runge-Kutta 2do Orden dy/dx

Referencia: Burden 5.4 p272 pdf282, Chapra 25.3 p740 pdf164, Rodriguez 9.1.7 p354, Boyce DiPrima 4Ed 8.4 p450

Como tema de introducción observar dos minutos del video sugerido a partir de donde se encuentra marcado el enlace. En este caso, en combate aereo, las armas se encuentran fijas en las alas.

Video Revisar:

Luego de observar el video de introducción conteste las siguientes preguntas:
¿ Que trayectoria siguen los proyectiles al salir del cañon?
¿ Que trayectoria siguen los aviones, el perseguido y el que caza?
¿ Cuál es la relación entre las trayectorias de los dos aviones?


Los métodos de Runge-Kutta  mejoran la aproximación a la respuesta sin requerir determinar las expresiones de las derivadas de orden superior. Los métodos usan una corrección a la derivada tomando valores de puntos alrededor referenciado al tamaño de paso h.

Por ejemplo, Runge-Kutta de 2do Orden usa el promedio entre los incrementos xi y xi+h, calculados como términos K1 y K2.

# EDO. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano   = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    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 + (K1+K2)/2
        xi = xi + h
        
        estimado[i] = [xi,yi]
    return(estimado)

Ejercicio

Para probar el algoritmo se usa la ecuación del problema presentado en  ‘EDO con Taylor‘ :

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

que aplicado con Runge Kutta, se obtiene:

estimado[xi,yi]
[[ 0.          1.        ]
 [ 0.1         1.2145    ]
 [ 0.2         1.4599725 ]
 [ 0.3         1.73756961]
 [ 0.4         2.04856442]
 [ 0.5         2.39436369]]
Error máximo estimado:  0.00435758459732
entre puntos: 
[ 0.          0.00067092  0.00143026  0.0022892   0.00326028  0.00435758]
>>> 

Compare los resultados con Taylor de 2 y 3 términos.

Los resultados se muestran también en la gráfica:

Se adjunta el programa de prueba que usa la función rungekutta2(d1y,x0,y0,h,muestras)  :

# PROGRAMA PRUEBA
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1

# INGRESO
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: y -x**2 + x + 1
x0 = 0
y0 = 1
h  = 0.1
muestras = 5

# PROCEDIMIENTO
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xi = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# SALIDA
print('estimado[xi,yi]')
print(puntosRK2)

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

yi_psol  = y_sol(xi)
errores  = yi_psol - yiRK2
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)

# Gráfica
import matplotlib.pyplot as plt

plt.plot(xis,yis, label='y conocida')
plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[x0,y0]')
plt.plot(xi[1:],yiRK2[1:],
         'o',color='m',
         label ='y Runge-Kutta 2 Orden')

plt.title('EDO: Solución con Runge-Kutta 2do Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

6.1 EDO con Taylor

Referencia: ejemplo Rodriguez 9.1.1 p335. Chapra 25.1.3 p731 pdf755

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 la derivada correspondiente.

Ejemplo:

Se pide 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

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

La solución empieza usando la Serie de Taylor para tres términos:

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)

Luego, al despejar de la fórmula planteada el valor de  y’, se puede obtener y" a combinar las ecuaciones.

y' = y -x^2 +x +1 y'' = y' -2x + 1 = (y -x^2 +x +1) -2x + 1 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.

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


Algoritmo con Python

Para simplificar los calculos 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 3 términos 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i,0:2] = [x,y]
    return(estimado)

La función se puede usar en un programa que plantee resolver el problema propuesto:

# PROGRAMA PRUEBA
# 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
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# SALIDA
print('estimado[xi, yi, d1yi, d2yi]')
print(puntos)

obteniendo como resultados:

estimado[xi,yi]
[[ 0.          1.        ]
 [ 0.1         1.215     ]
 [ 0.2         1.461025  ]
 [ 0.3         1.73923262]
 [ 0.4         2.05090205]
 [ 0.5         2.39744677]]

Cálculo de Error

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

# 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)

con los siguientes resultados:

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

Gráfica

Para visualizar los resultados se usan los puntos estimados y muchos más puntos para la solución conocida en el intervalo de estudio. Con la gráfica se podrán observar las diferencias entre las soluciones encontradas.

# 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)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xis,yis, label='y conocida')
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()

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]

5.4.1 Cuadratura con dos puntos – Experimento

Para visualizar el concepto de Cuadratura de Gauss de dos puntos considere lo siguiente:

Se tiene un corte transversal del un recipiente rectangular lleno de líquido limitado en x entre [-1,1], al que se le ha depositado encima otro recipiente con perfil f(x) = x2 hasta que reposan sus extremos en x[-1,1].

La altura de ambos recipientes es la misma.

La superficie entre f(x) y el eje x es el integral de f(x) en el intervalo.

Si suponemos que la figura es el corte transversal de una vasija y la parte en amarillo es líquida, la vasija ha desplazado el líquido que ocupa ahora el «área» mostrada en la gráfica que corresponde al integral de f(x)=x2. entre [-1,1].

Ahora, suponga que se perfora el perfil de f(x) en dos puntos equidistantes cercanos  x=0. Los orificios permitirían el desplazamiento del liquido al interior de f(x) que dejando pasar suficiente tiempo, perimitiría tener todo el líquido en el recipiente rectangular entre [-1,1] como una línea horizontal.

Podría medir la altura que tiene el líquido y que tiene un equivalente en un punto f(x1). Debería encontrar el valor de x1 que permite disponer del mismo valor entre el área debajo de f(x) y el rectángulo del corte transversal amarillo ahora formado.

Se usa el resultado analítico del integral restando el área del rectángulo obtenido al evaluar la funcion f(x) entre [0,1], teniendo un problema de busqueda de raíces. Obtenemos el valor de x1.

Se muestra que el área bajo f(x) es equivalente al área del rectángulo conformado.

Si  utilizamos el desplazamiento horizontal desde el centro para un punto encontrado como un «factor», tendremos que el área del rectángulo se mantendría equivalente, y el desplazamiento proporcional a la mitad del intervalo si se varía el intervalo de observacion. Este factor coincide con el factor de Cuadratura de Gauss de dos puntos.

funcion  fx:   x**2
Integral Fx:   x**3/3
I analitico:   0.6666666666666666
I aproximado:  0.6666666666654123
desplaza centro:   0.5773502691890826
factor desplaza:   0.5773502691890826
Factor CuadGauss:  0.5773502691896258
erradoFactor:    1.2545520178264269e-12
error integral:  1.2543299732215019e-12

El error del integral es del orden de 10-12


Cambiamos la figura geométrica a un trapecio generado por la recta que pasa por los puntos xi desplazados desde el centro.

Usamos la función f(x) = x2 + x + 1, observaremos si los resultado son equivalentes.

La figura al inicio del experimento será:

Luego de realizar realizar el mismo cálculo anterior usando un equivalente a trapecio se tiene:

con valores numéricos:

funcion  fx:   x**2 + x + 1
Integral Fx:   x**3/3 + x**2/2 + x
I analitico:   2.6666666666666665
I aproximado:  2.6666666666654124
desplaza centro:   0.5773502691890826
factor desplaza:   0.5773502691890826
Factor CuadGauss:  0.5773502691896258
erradoFactor:    1.2545520178264269e-12
error integral:  1.2541079286165768e-12

El error del integral es también del orden de 10-12, además observe que el factor de cuadratura de Gauss se mantiene.


Tarea

Realice el experimento usando un polinomio de grado superior y observe los errores para el integral y las diferencia con el coeficiente de Cuatratura de 2 puntos.


Algoritmo con Python

Para resumir la cantidad de instrucciones, se usa el método de la bisección desde la librería scipy y el subgrupo de funciones de optimización.

Los cálculos para realizar las gráficas se tratan en un bloque luego de mostrar los resultado principales.

# Integración: Cuadratura de Gauss de dos puntos
# modelo con varios tramos entre [a,b]
# para un solo segmento.
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import scipy.optimize as op

# INGRESO
x = sym.Symbol('x')

# fx = (x)**2
fx = x**2 + x + 1
# fx = 0.2 + 25.0*x-200*(x**2)+675.0*(x**3)-900.0*(x**4)+400.0*(x**5)

a = -1
b = 1
muestras = 51
tolera = 1e-12
iteramax = 100

# PROCEDIMIENTO
# Desarrollo analítico con Sympy
Fx = sym.integrate(fx,x)
Fxn = sym.lambdify('x',Fx,'numpy')
Fiab = Fxn(b)-Fxn(a)

# Busca igualar trapecio con Integral analitico
fxn = sym.lambdify('x',fx,'numpy')
base = b-a
mitad = base/2
xc = (a+b)/2  # centro
diferencia = lambda x: Fiab-base*(fxn(xc-x)+fxn(xc+x))/2
desplazado =  op.bisect(diferencia,0,mitad,
                        xtol=tolera,maxiter=iteramax)
factor = desplazado/mitad

# Integral aproximando con trapecio
x0 = xc - factor*mitad
x1 = xc + factor*mitad
Faprox = base*(fxn(x0)+fxn(x1))/2

# Integral cuadratura Gauss
xa = xc + mitad/np.sqrt(3)
xb = xc - mitad/np.sqrt(3)
FcuadG = base*(fxn(xa)+fxn(xb))/2
erradofactor = np.abs(FcuadG - Faprox)
erradoIntegral = np.abs(Fiab-Faprox)
# SALIDA
print('funcion  fx:  ', fx)
print('Integral Fx:  ', Fx)
print('I analitico:  ', Fiab)
print('I aproximado: ', Faprox)
print('desplaza centro:  ', desplazado)
print('factor desplaza:  ', factor)
print('Factor CuadGauss: ', 1/np.sqrt(3))
print('erradoFactor:   ', erradofactor)
print('error integral: ', erradoIntegral) 

# Grafica
# Para GRAFICAR 
# Para gráfica f(x)
xi = np.linspace(a,b,muestras)
fi = fxn(xi)

# Para gráfica Trapecio
m = (fxn(x1)-fxn(x0))/(x1-x0)
trapeciof = lambda x: fxn(x0)+m*(x-x0)
trapecioi = trapeciof(xi)

# Areas Trapecio para cada punto que busca
k = int(muestras/2)
xicg = xi[k:muestras-1]
Fcg = [base*(fxn(xi[k+0])+fxn(xi[k-0]))/2]
for i in range(1,k,1):
    untrapecio = base*(fxn(xi[k+i])+fxn(xi[k-i]))/2
    Fcg.append(untrapecio)

# Punto buscado
Fiaprox = base*(fxn(x1)+fxn(x0))/2

Fi = Fxn(xi)-Fxn(a)

# Areas de curvas y trapecio

plt.subplot(211) # Grafica superior
plt.xlim(a,b)
plt.plot(xi, fi, label='f(x)')
# Solo fi
# plt.fill_between(xi,0, fi,
#                label='integral fi',
#                 color='yellow')
# usando cuadratura
plt.fill_between(xi,0, trapecioi,
                 label='Cuadratura 2 puntos',
                 color='yellow')
plt.axvline(x0,color='white')
plt.axvline(x1,color='white')
plt.plot([x0,x1],[fxn(x0),fxn(x1)],
         'ro', label='x0,x1')
plt.axvline(0,color='black')
plt.xlabel('x')
plt.ylabel('f(x) y Cuadratura de 2 puntos')
plt.legend()

# Valores de integrales
plt.subplot(212) # Grafica inferior
plt.xlim(a,b)
plt.axhline(Fiab, label='F[a,b]')
# plt.plot(xi,Fi,label='F(x)')
plt.plot(xicg,Fcg,color='orange',label='Aprox Trapecio')

plt.axvline(x1,color='yellow')
plt.axvline((1/np.sqrt(3))*(b-a)/2 + xc ,color='magenta')
plt.plot(x1,Fiaprox,'ro', label='x0,x1')
plt.axvline(0,color='black')

plt.xlabel('x')
plt.legend()
plt.ylabel('Integrando')
plt.show()

5.4 Cuadratura de Gauss con Python

Referencia: Chapra 22.3 p655 pdf679, Rodriguez 7.3 p294, Burden 4.7 p220 pdf230

La cuadratura de Gauss aproxima el integral de una función en un intervalo [a,b] centrado en cero mediante un cálculo numérico con menos operaciones y evaluaciones de la función. Se representa como una suma ponderada:

I \cong c_0f(x_0) + c_1f(x_1)

para la fórmula de dos puntos se tiene obtiene:

c_0 = c_1 = 1 x_0 = -\frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}

para un intervalo de evaluación desplazado en el eje x se requiere convertir los puntos al nuevo rango. Se desplaza el punto cero al centro del intervalo [a,b] y se obtiene:

x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{b-a}{2}(f(x_a) + f(x_b))

cuya fórmula es semejante a una mejor aproximación de un trapecio, cuyos promedios de alturas son puntos internos de [a,b], concepto mostrado en la gráfica.


Ejercicio

El algoritmo se desarrolla en un tramo en el intervalo [a,b] junto a la gráfica para mostrar el concepto. Para el ejemplo el integral buscado es:

\int_1^3 \sqrt{x} \sin(x) dx

que usando cuadratura de Gauss con un solo intervalo[a,b] tiene como resultado:

Integral:  2.03821975687
>>>

el resultado se puede mejorar aumentando el número de tramos en el intervalo [a,b]. Por ejemplo, el resultado usando 4 tramos el resultado es semejante al usar el método del trapecio con 128 tramos, lo que muestra el ahorro en calculos entre los métodos

Integral:  2.05357719003
>>> 

Algoritmo con Python

Instrucciones usando la cuadratura de Gauss como una función

# Integración: Cuadratura de Gauss de dos puntos
# modelo con varios tramos entre [a,b]
import numpy as np
import matplotlib.pyplot as plt

# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    return(area)

# INGRESO
fx = lambda x: np.sqrt(x)*np.sin(x)

# intervalo de integración
a = 1
b = 3
tramos = 4

# PROCEDIMIENTO
muestras = tramos+1
xi = np.linspace(a,b,muestras)
area = 0
for i in range(0,muestras-1,1):
    deltaA = integraCuadGauss2p(fx,xi[i],xi[i+1])
    area = area + deltaA
# SALIDA
print('Integral: ', area)

Gráfica por tramos

La gráfica que complementa el resultado anterior, se realiza añadiendo las intrucciones presentadas a continuación.

Considere que la gráfica es útil con pocos tramos en el intervalo[a,b]

# GRAFICAR por cada Segmento
# para concepto con pocos segmentos
x0 = -1/np.sqrt(3) 
x1 = 1/np.sqrt(3)

# arregos para gráficas
xi = np.array([])
fi = np.array([])

xat = np.array([])
xbt = np.array([])

recta = np.array([])

muestrastramo = 10
subtramo = np.linspace(a,b,muestras)

for i in range(0,tramos,1):
    at = subtramo[i]
    bt = subtramo[i+1]
    
    xit = np.linspace(at,bt,muestrastramo)
    fit = fx(xit)
    
    xi = np.concatenate((xi,xit))
    fi = np.concatenate((fi,fit))

    # puntos xa y xb por tramo
    xa = (bt+at)/2 + (bt-at)/2*(x0)
    xb = (bt+at)/2 + (bt-at)/2*(x1)
    
    xat = np.concatenate((xat,[xa]))
    xbt = np.concatenate((xbt,[xb]))
    
    # Recta entre puntos x0 y x1 por tramo
    m = (fx(xb)-fx(xa))/(xb-xa)
    b0 = fx(xa)- m*xa
    linea = b0 + m*xit
    recta = np.concatenate((recta,linea))

# Marcadores 'o' de xa y xb por tramos
puntox = np.concatenate((xat,xbt))
puntoy = fx(puntox)

# Graficando
# Trazado de lineas
plt.plot(xi,recta, label = 'grado 1', color = 'tab:orange')
plt.fill_between(xi,0,recta, color='tab:olive')
plt.plot(xi,fi, label='f(x)', color = 'blue')

# Verticales para dividir los tramos
for j in range(0,len(subtramo),1):
    plt.axvline(subtramo[j], color='tab:gray')
    
# Marcadores de puntos xa y xb por tramos
for j in range(0,len(xat),1):
    plt.axvline(xat[j], color='w')
    plt.axvline(xbt[j], color='w')
    
plt.plot(puntox,puntoy, 'o', color='g')

plt.title('Integral: Cuadratura Gauss')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()

plt.show()