Sistemas de ecuaciones diferenciales ordinarias

Estudiamos sistemas de ecuaciones diferenciales no lineales: $$ y = f(t,y) $$ donde $y$ es una función vectorial (toma valores en $\mathbb{R}^d$) de una sóla variable real $t$, que se suele identificar con el tiempo, pero también puede ser una variable espacial (por ejemplo si estudiamos la tensión a lo largo de una viga). Para hacer énfasis es que $y$ es una función del tiempo, podemos escribir $$ y(t) = f(t,y(t)) $$

Ecuación funcional

En la ecuación anterior $y = f(t,y)$, la incógnita es una función $y$. No sólo no conocemos el valor de la función $y(t)$ para cada instante temporal $t$, sino que además no sabemos para qué valores de $t$ está definida (en otras palabras, no conocemos su dominio).

Problemas de valor inicial

En este cuaderno, estudiamos problemas de valor inicial (PVI), donde la variable $t$ de la que depende la función $y(t)$ es el tiempo, y conocemos el valor de $y$ en el instante inicial $t_0$: $$\begin{array}{lll} y (t) & = & f (t, y (t))\\ y (t_0) & = & y_0 \end{array} $$

Método de Euler

Habéis visto en clase el método de Euler, que aproxima la solución al PVI anterior mediante la iteración: $$\begin{array}{lll} t_k &=& t_0 + h\cdot k\\ y_{k+1} & = & y_k + h\cdot f (t_k, y_k)\\ y_0 & = & y(t_0)=y_0 \end{array}$$ Este método converge, y el error de truncamiento al aproximar $y(t)$ usando un paso $h$ es: $$ |{\text{Error}_{\text{trunc}}}|\leq {\frac {hM}{2L}}(e^{L(t-t_{0})}-1)\qquad \qquad $$ donde $L$ es la constante Lipschitz de $f$ en su segundo argumento. Si $f$ es de clase $C^1$, entonces $L$ es el máximo de la derivada parcial $\frac{\partial f}{\partial y}$.

In [9]:
function [ts, ys] = Euler(f,y0,a,b,n)

    % La f de entrada debe ser introducida en la forma @(t,y)f(t,y)

    % La EDO a resolver será de la forma y'=f(t,y)
    % (a,b) será el intervalo
    % y(a)=alpha, será la condición inicial del PVI
    % n será el número de iteraciones o de particiones del intervalo

    %%%%%%%%% IMPLEMENTACIÓN DEL MÉTODO DE EULER %%%%%%%%%%%%%%%%%%%
    %% PASO 1: definimos el estado inicial de los vectores "st" y "sw" 
    %% donde se guardarán los valores de t y w respectivamente.
    ts=[a];
    ys=[y0'];
    %% PASO 2: definimos el tamaño de paso h
    h = (b - a)/n;
    %% PASO 3: inicializamos los valores de las variables de iteración  t y w 
    %% en los datos iniciales
    t=a; 
    w=y0;
    %% PASO 4: introducimos la ley de recurrencia de Euler
    for i=1:n
        w=w+h*f(t,w);
        t=t+h;
       %% cargamos los nuevos valores de t y w obtenidos en los vectores 
       %% "st" y "sw", respectivamente
       ts=[ts t];
       ys=[ys w'];
    end
end

Ejemplo: péndulo amortiguado

Un péndulo amortiguado sigue la ecuación:

$$ m \cdot l\cdot \theta''+k \cdot l\cdot \theta'+m \cdot g \cdot \sin(\theta)=0 $$

donde

  • $l$ longitud del péndulo (= 10 m).
  • $\theta$ ángulo que forma la cuerda con la vertical.
  • g aceleración de la gravedad (≈ 10 m/sg ).
  • m masa de la bola (= 20 kg).
  • k coeficiente de fricción del medio (= 4 kg/sg).

El péndulo comienza en la posición de reposo. Se imprime al péndulo una velocidad inicial de 1 rad/s. Calculamos la trayectoria del péndulo mediante del método de Euler.

Para poder resolverla con el método de Euler, necesitamos escribir esta ecuación de segundo orden como un sistema de ecuaciones diferenciales de primer orden, introduciendo una nueva variable $v$: $$\left\{\begin{array}{lll} \theta' & = & v\\ v' & = & - \frac{k}{m} v - \frac{g}{l} \sin (\theta) \end{array}\right. $$ con las condiciones iniciales $$\left\{\begin{array}{lll} \theta(t_0) & = & 0\\ v(t_0) & = & 1 \end{array}\right.$$

In [12]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2), -(k/m)*y(2)-(g/l)*sin(y(1))]);

t0 = 0;
tf = 20;
# y = (theta, v)
y0 = [0,1];
N = 120;

[ts, ys] = Euler(fun, y0, t0, tf, N);
thetas = ys(1,:);
vs = ys(2,:);

plot(ts, thetas)
xlabel('t')
ylabel('theta')
title('Las oscilaciones se amortiguan...')

Curiosidad: observa qué ocurre al aumentar la velocidad inicial.

In [22]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2), -(k/m)*y(2)-(g/l)*sin(y(1))]);

t0 = 0;
tf = 20;
# y = (theta, v)
y0 = [0,3];
N = 120;

[ts, ys] = Euler(fun, y0, t0, tf, N);
thetas = ys(1,:);
vs = ys(2,:);

plot(ts, thetas)
xlabel('t')
ylabel('theta')
title('El impulso inicial hace dar al pendulo una vuelta entera, \n y luego las oscilaciones se amortiguan...')

Usando un método "profesional"

Usamos el método ode45, que usa un método muy popular: el método adaptativo Runge-Kutta 45:

  • Combina un método de Runge-Kutta de orden $4$ con otro de orden $5$ para estimar el error.
  • Si el error es demasiado grande, reduce el paso $h$. Si es muy pequeño, lo aumenta.
  • Sigue avanzando con este método de Runge-Kutta de orden $5$ para el paso elegido.

Hablaremos de métodos adaptativos más despacio en otro cuaderno, por ahora nos contentaremos con observar una dificultad al usar estos métodos: El paso es de longitud variable y a menudo los pasos son grandes, porque el método no necesita un paso pequeño para obtener un error reducido.. que nos permite elegir uno entre varios métodos numéricos.

# Intervalo de tiempo
t_span = [0,20]
# Punto inicial
y0 = [0,1]

sol = integ.solve_ivp(
    fun, t_span, y0, 
    method='RK45')

Aunque ode45 es adecuado para muchos problemas, no siempre es el mejor método. Consulta la documentación de matlab sobre cómo elegir un método numérico.

In [23]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2), -(k/m)*y(2)-(g/l)*sin(y(1))]);

tspan = [0, 20];
# y = (theta, v)
y0 = [0,3];
N = 120;

[ts, ys] = ode45(fun, tspan, y0);
thetas = ys(:,1);
vs = ys(:,2);

plot(ts, thetas)
xlabel('t')
ylabel('theta')
title('El impulso inicial hace dar al pendulo una vuelta entera, \n y luego las oscilaciones se amortiguan...')

El método ode45 parece tener más error que el método de Euler, pero es una ilusión. Lo que ocurre es que los instantes de tiempo que necesita para encontrar la solución aproximada están muy separados.

Si queremos una trayectoria suave, es necesario interpolar. Podríamos hacerlo a mano usando alguna de las técnicas que hemos aprendido, pero ode45 trae de serie un algoritmo de interpolación que está especialmente adaptado al método numérico elegido. Para ode45, como tenemos aproximaciones a $f$ y a todas sus derivadas hasta orden 4, solve_ivp utiliza interpolación polinómica de Hermite a trozos con polinomios de grado 4 a trozos.

La forma de pedir a ode45 que interpole es pasarle una lista de tiempos en vez de sólo el intervalo tspan.

In [27]:
m = 20;
l = 10;
k = 4;
g = 10;

fun = @(t,y)([y(2), -(k/m)*y(2)-(g/l)*sin(y(1))]);

N = 200;
t0 = 0;
tf = 20;
tspan = linspace(t0, tf,N);
# y = (theta, v)
y0 = [0,3];
N = 120;

[ts_ode45, ys_ode45] = ode45(fun, [t0, tf], y0);

[ts, ys] = ode45(fun, tspan, y0);
thetas = ys(:,1);
vs = ys(:,2);

hold on;
plot(ts_ode45, ys_ode45(:,1), 'ro')
plot(ts, thetas, 'b-')
hold off;
xlabel('t')
ylabel('theta')
title('El impulso inicial hace dar al pendulo una vuelta entera, \n y luego las oscilaciones se amortiguan...')

Ejercicio

  • En vez de dibujar la posición con respecto al tiempo, dibuje la evolución de la energía cinética respecto al tiempo.
  • Dibuja el diagrama de fases del péndulo no lineal. Compara el resultado con el diagrama de fases del péndulo lineal que pedimos en el cuaderno de NIVELACION/edos_lineales.
In [ ]:
In [ ]:
In [ ]:

Ejercicio

La ecuación de van der Pol modeliza ciertos circuitos eléctricos, entre otras cosas:

The Van der Pol oscillator was originally proposed by the Dutch electrical engineer and physicist Balthasar van der Pol while he was working at Philips. Van der Pol found stable oscillations, which he subsequently called relaxation-oscillations and are now known as a type of limit cycle in electrical circuits employing vacuum tubes.

La ecuación para la corriente que fluye por cierta parte del circuito es: $$ {\displaystyle {d^{2}x \over dt^{2}}-\mu (1-x^{2}){dx \over dt}+x=0.} $$ Se pide:

  • Para $\mu=0.5$, resolver el problema de valor inicial $x(0)=1, x'(0)=0$ con el método de Euler con h=0.1 y con algún método adaptativo, hasta alcanzar $t=2$. Comparar la precisión y el "número de evaluaciones de la ecuación f", para ambos métodos.
  • Realice un dibujo de plano de fases para la solución del sistema, integrando hasta $t=10$ con un método de su elección.
In [ ]:
In [ ]:
In [ ]: