Métodos numéricos de integración de EDOs

from pylab import *

Consideremos el pvi $$ \left\{\begin{aligned} &\dot{x}=1+t-2x\\ &x(0)=1, \end{aligned}\right. $$ que tiene por solución $$ x=\frac{1}{4}(1+2t+3e^{-2t}). $$

# Ecuación diferencial x' = f(t,x)
def f(t, x):
    return 1 + t - 2*x

# Solución exacta
t = linspace(0, 1);
x = (1 + 2*t + 3*exp(-2*t)) / 4.0
plot(t,x,'g');

Método de Euler (progresivo)

Para simplificar tomo todos los $t_i$ a la misma distancia: $t_i-t_{i-1}=h$.

A modo de calculadora

h = 0.2

t0 = 0
x0 = 1

x1 = x0 + h*f(t0,x0)
t1 = t0+h

x2 = x1 + h*f(t1,x1)
t2 = t1+h

x3 = x2 + h*f(t2,x2)
t3 = t2+h

x4 = x3 + h*f(t3,x3)
t4 = t3+h

x5 = x4 + h*f(t4,x4)
t5 = t4+h

tiempos = [t0,t1,t2,t3,t4,t5]
X = [x0,x1,x2,x3,x4,x5]

plot(tiempos, X, 'bo:', t, x, 'g');

Definamos una función que implemente el método de Euler. Queremos usarla en la forma

Euler_p(fun, t_span=[0,1], x0=1, N=5) ---> tiempos, X

# Método Euler progresivo
def Euler_p(fun, t_span, x0, N):
    (a, b) = t_span
    x = asarray(x0)
    tiempos, h = linspace(a, b, N+1, retstep=True)
    X = [x]
    for t in tiempos[0:-1]:
        x = x + h*fun(t,x)
        X.append(x)
    return (tiempos, asarray(X).T)

r, X = Euler_p(fun=f, t_span=[0,1], x0=1, N=5)
plot(r, X, 'bo:',  t, x, 'g');

r , X = Euler_p(fun=f, t_span=[0,1], x0=1, N=50)
plot(r, X, 'b',  t, x, 'g');

Otros ejemplos

La ecuación logística $\dot{x}=rx(1-x)$.

def logistica(t,x):
    r=0.3
    return r*x*(1-x)

r,X = Euler_p(fun=logistica, t_span=[0,100], x0=0.00001, N=50)
plot(r, X, 'b');

La ecuación $\dot{x}=x^2$

def xdos(t,x):
        return x**2
r, X = Euler_p(fun=xdos, t_span=[0,1], x0=1, N=500)
plot(r, X, 'b');

Sistemas de ecuaciones

El programa anterior vale tanto para ecucaiones escalares como para sistemas. Veamos como resolver el pvi $$ \left\{\begin{aligned} &\dot{x}=-3x+2y\\ &\dot{y}=-2x-3y\\ &x(0)=0\\ &y(0)=1 \end{aligned}\right. $$ para valores de $t$ en el intervalo $[0,2]$

def f(t,u):
    x,y = u
    return array([-3*x+2*y,-2*x-3*y])

t, U = Euler_p(fun=f, t_span=[0,2], x0=[0,1], N=8)
# La columna i-ésima de U es [x_i,y_i]
# La primera fila de U son las x, la segunda son las y.
U

array([[ 0.        ,  0.5       ,  0.25      , -0.03125   , -0.09375   ,
        -0.03710938,  0.01074219,  0.01696777,  0.00512695],
       [ 1.        ,  0.25      , -0.1875    , -0.171875  , -0.02734375,
         0.04003906,  0.02856445,  0.00177002, -0.00804138]])

t,U = Euler_p(fun=f, t_span=[0,2],x0=[0,1],N=100)
x,y = U

figure()
xlabel('t')
plot(t,x, label='x');
plot(t,y, label='y');
legend()

figure()
xlabel('x')
ylabel('y')
plot(x,y,'r');

Usando las librerías de Python

Python tiene ya programados multitud de métodos, muy eficientes, para la integración numérica de EDOs. La manera más cómoda de usarlos en mediante la función solve_ivp.

La notación usada por solve_ivp es $$ \left\{\begin{align}&y'=f(t,y)\\&y(t_0)=y_0,\end{align}\right. $$ es decir, nuestra variable $x$ él la llama $y$.

Devuelve un objeto python con mucha información (que nosotros no necesitamos).

Veamos como usarlo.

from scipy.integrate import solve_ivp

t = linspace(0,2,9)
result = solve_ivp(fun=f, t_span=[0,2], y0=[0,1], t_eval=t)
#result

# result.y es la matriz con los valores de la variable dependiente
result.y

array([[ 0.        ,  0.22650781,  0.18778725,  0.10509696,  0.04521333,
         0.01403154,  0.00154069, -0.00185156, -0.00187785],
       [ 1.        ,  0.41451869,  0.12045908,  0.00735179, -0.0207666 ,
        -0.01885086, -0.01098926, -0.00490266, -0.00161167]])

t = linspace(0,2,100)
result = solve_ivp(fun=f, t_span=[0,2], y0=[0,1], t_eval=t)
x,y = result.y

figure()
xlabel('t')
plot(t,x, label='x');
plot(t,y, label='y');
legend()

figure()
xlabel('x')
ylabel('y')
plot(x,y,'r');

Métodos gráficos

Dibujo del campo vectorial y diagrama de fases

Se utiliza la función quiver (en castellano carcaj, bolsa para las flechas)

def N(x,y):
    return sqrt(x**2+y**2)

def F(x,y):
    return array([-y,x])/N(x,y)

X,Y = meshgrid( linspace(-1,1,10), linspace(-2,2,20) )
U,V = F(X,Y)
color = N(X,Y)
quiver(X,Y,U,V,color,pivot='mid');

streamplot(X,Y,U,V, density=0.7);

Utilidad para no repetir código

def plot_vector_field(f, x_range, y_range, N=[30,30]):
    a,b = x_range
    c,d = y_range
    nx,ny = N
    X,Y = meshgrid( linspace(a,b,nx), linspace(c,d,ny) )
    U,V = f([X,Y])
    return quiver(X,Y,U,V, pivot='mid')

def F(u):
    x,y = u
    dx = 3*x+2*y
    dy = (2*y-x)*x
    return array([dx,dy])/sqrt(dx**2+dy**2)

plot_vector_field(F,[-1,1],[-1,1],[20,30]);

Conjuntos de nivel de una integral primera

Para la ecuación diferencial $$ (t+x)\,dx+x\,dt $$ la función $$ F(t,x)= \frac{1}{2}x^2+tx $$ es la función potencial, una integral primera. Las soluciones satisfacen $F(t,x(t))=c$, es decir, están contenidas en los conjuntos de nivel de $F$.

La gráfica de los conjuntos de nivel se obtiene con contour o contourf.

def F(t,x):
    return x**2+2*t*x

T, X = meshgrid( linspace(-1.5, 1.5,101), linspace(-1.5, 1.5,101) )
Z=F(T,X)

contour(T, X, Z, 20);

contourf(T, X, Z, 20);