Estas ecuaciones se pueden escribir vectorialmente como
$$ \frac{d \vec{y}}{d\xi} = \vec{f}\left(\vec{\xi},\,\,\vec{y}\right) $$$$ \vec{y} = \left(\begin{matrix} g \\ y_2 \end{matrix}\right),\,\, \vec{f} = \left(\begin{matrix} y_2 \\ T \end{matrix}\right) $$defs:
El problema con este algoritmo es que es inestable, y por lo tanto los errores asociados crecen en vez de cancelarse. Ahora vamos a analizar la estabilidad del algoritmo.
Sea $\bar{y}$ la solucion exacta al algoritmo numerico ($\neq$ a la solucion del problema)
$$\bar{y}_{n+1} = \bar{y}_{n-1} + 2hf\left(x_n,\,\,\bar{y}_n\right) \qquad\qquad (**)$$Sea $\epsilon_n$ el error de maquina asociado al paso $n$
$$\implies y_n = \bar{y}_n + \epsilon_n$$Esto se reemplaza en $(*)$ y se obtiene: $$\bar{y}_{n+a}+\epsilon_{n+1} = \bar{y}_{n-1}+\epsilon_{n-1}+2f\left(x_n,\bar{y}_n+\epsilon_n\right)h$$
Haciendo una expansion de taylor de $f$ en torno a $\epsilon_n$:
$$f\left(x_n,\,\,\bar{y}_n+\epsilon_n\right)\sim f\left(x_n,\,\,\bar{y}_n\right) + \frac{\partial f}{\partial y}\biggr\vert_{\bar{y}_n}\epsilon_n$$Esto despues se reemplaza en la expresion original, y aprovechando que la solucion exacta del algoritmo es $(**)$ la ecuacion queda:
$$\epsilon_{n+1} = \epsilon_{n-1} + 2h\frac{\partial f}{\partial y}\biggr\vert_{\bar{y}_n}\epsilon_n$$El ultimo termino de esta expresion se puede asumir constante y se renombra $2h\frac{\partial f}{\partial y} = 2\gamma = $ const
Suponiendo que $\epsilon_n = \epsilon_0\lambda^n$
$$\epsilon_0\lambda^{n+1} = \epsilon_0\lambda^{n-1} + 2\gamma\epsilon_0\lambda^n$$$$\lambda^2 - 2\gamma\lambda - 1 = 0$$$$\lambda = \gamma\pm\sqrt{1+\gamma^2}$$La solucion general es:
$$\epsilon_n = A\lambda_{+}^n + B\lambda_{-}^n$$Se tienen los siguientes casos: $$\gamma > 0 \implies \lambda_{+} > 1 \implies \epsilon_n \rightarrow \infty\quad \text{as}\quad n \rightarrow \infty$$
$$\gamma < 0 \implies \lambda_{-} < 1- \implies \epsilon_n \rightarrow \pm\infty\quad \text{as}\quad n \rightarrow \infty$$Este algoritmo suele ser utilizado para calcular derivadas de primer orden, o resolver ecuaciones en donde el error asociado al problema es menor al error acumulado de maquina.
Buscamos un algoritmo mas estable...
Reemplazando 2) en 1) $$y_{n+1} = y_n + hy'_{n+1/2} + O\left(h^3\right)$$
Con $f_n = y'_n$
Con $\phi_0$ algun angulo, $\omega_0 = 0 \implies \phi(t) = A\cos(\lambda t + \delta) \quad \dot{\phi} = -A\lambda\sin(\lambda t + \delta)\quad$ con $\lambda = \sqrt{\frac{g}{l}}$
In [ ]:
"""
Implementa la solucion de RK2 para el problema de un pendulo de largo 1
"""
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# Constantes
g = 9.8
# Condiciones iniciales
phi_t0 = np.pi / 90.
w_t0 = 0.
# Solucion para pequenas oscilaciones
freq = np.sqrt(g)
t = np.linspace(0, 4 * 2 * np.pi / freq, 10000) # 4 ciclos
phi_po = phi_t0 * np.cos(freq * t)
plt.plot(t, phi_po, label='pequenas oscilaciones')
plt.legend(loc=8)
Out[2]:
In [3]:
# Usando RK2
# h es el tamano del paso para cada iteracion!!
def f(phi, omega):
output = [omega, -g * np.sin(phi)]
return output
def calc_k1(f, phi_n, omega_n, h):
fn = f(phi_n, omega_n)
output = [h * fn[0], h * fn[1]]
return output
def calc_k2(f, phi_n, omega_n, h):
k1 = calc_k1(f, phi_n, omega_n, h)
f_mid = f(phi_n + k1[0]/2, omega_n + k1[1]/2)
output = [h * f_mid[0], h * f_mid[1]]
return output
def rk2_step(f, phi_n, omega_n, h):
k2 = calc_k2(f, phi_n, omega_n, h)
phi_nxt = phi_n + k2[0]
omega_nxt = omega_n + k2[1]
output = [phi_nxt, omega_nxt]
return output
Nsteps = 10000
h = 4 * 2 * np.pi / freq / Nsteps # tiempo / pasos
phi_arr = np.zeros(Nsteps)
omega_arr = np.zeros(Nsteps)
# Condiciones iniciales
phi_arr[0] = phi_t0
omega_arr[0] = w_t0
for i in range(1, Nsteps):
phi_arr[i], omega_arr[i] = rk2_step(f, phi_arr[i-1], omega_arr[i-1], h)
plt.plot(t, phi_arr, label='rk2')
plt.legend(loc=8)
Out[3]:
In [4]:
# Comparacion de ambos resultados
plt.plot(t, phi_po, label='pequenas oscilaciones')
plt.plot(t, phi_arr, label='rk2', color='green')
plt.legend(loc=8)
Out[4]:
El algoritmo de Runge-Kutta 4 es mas preciso que los anteriores, pero, es estable? La repuesta a esa pregunta depende fuertemente del problema a resolver.
A continuacion se analiza la estabilidad de Runge-Kutta 4 para un problema determinado.
Consideremos $$\frac{dx}{dy} = \lambda y$$ Supongamos que $$ y_n = \bar{y}_n + \epsilon_n$$
Luego
$$ k_1 = hf\left(x_n,\,\, y_n\right) = h\lambda y_n $$Para que el algoritmo sea estable se debe cumplir que $-2.7853 < h\lambda < 0$
Solucion analitica:
$$ x(t) = \sin(\sqrt{k}t) $$