Ecuaciones diferenciales

Abrir script en vivo

En este ejemplo se muestra cómo utilizar MATLAB® para formular y resolver varios tipos distintos de ecuaciones diferenciales. MATLAB ofrece varios algoritmos numéricos para resolver una gran variedad de ecuaciones diferenciales:

Problemas de valores iniciales
Problemas de valores de límites
Ecuaciones diferenciales con retardo
Ecuaciones diferenciales parciales

Problemas de valores iniciales

vanderpoldemo es una función que define la ecuación de van der Pol

$\frac{d^{2} y}{d t^{2}} - μ (1 - y^{2}) \frac{d y}{d t} + y = 0 .$

type vanderpoldemo

function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

La ecuación se escribe como un sistema de dos ecuaciones diferenciales ordinarias (ODE) de primer orden. Estas ecuaciones se evalúan para diferentes valores del parámetro $μ$ . Para una integración más rápida, debe elegir un solver apropiado basado en el valor de $μ$ .

Para $μ = 1$ , cualquiera de los solvers de ODE de MATLAB puede resolver eficazmente la ecuación de van der Pol. El solver ode45 es uno de esos ejemplos. La ecuación se resuelve en el dominio $[0, 20]$ con las condiciones iniciales $y (0) = 2$ y $\frac{dy}{dt} |_{t = 0} = 0$ .

tspan = [0 20];
y0 = [2; 0];
Mu = 1;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);

% Plot solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')

Figure contains an axes object. The axes object with title van der Pol Equation, mu blank = blank 1, xlabel t, ylabel solution y contains an object of type line.

Para magnitudes mayores de $μ$ , el problema se convierte en rígido. Esta etiqueta se aplica a problemas que resisten los intentos de ser evaluados con técnicas ordinarias. En su lugar, se necesitan métodos numéricos especiales para una integración rápida. Las funciones ode15s, ode23s, ode23t y ode23tb pueden resolver eficazmente problemas rígidos.

Esta solución de la ecuación de van der Pol para $μ = 1000$ utiliza ode15s con las mismas condiciones iniciales. Se necesita ampliar drásticamente el periodo de tiempo a $[0, 3000]$ para poder ver el movimiento periódico de la solución.

tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode15s(ode, tspan, y0);

plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000')
axis([0 3000 -3 3])
xlabel('t')
ylabel('solution y')

Figure contains an axes object. The axes object with title van der Pol Equation, mu blank = 1000, xlabel t, ylabel solution y contains an object of type line.

Problemas de valores de límites

bvp4c y bvp5c resuelven problemas de valores de límites para ecuaciones diferenciales ordinarias.

La función de ejemplo twoode tiene una ecuación diferencial escrita como un sistema de dos ODE de primer orden. La ecuación diferencial es

$\frac{d^{2} y}{d t^{2}} + | y | = 0 .$

type twoode

function dydx = twoode(x,y)
%TWOODE  Evaluate the differential equations for TWOBVP.
%
%   See also TWOBC, TWOBVP.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

dydx = [ y(2); -abs(y(1)) ];

La función twobc tiene las condiciones de límites para el problema: $y (0) = 0$ e $y (4) = - 2$ .

type twobc

function res = twobc(ya,yb)
%TWOBC  Evaluate the residual in the boundary conditions for TWOBVP.
%
%   See also TWOODE, TWOBVP.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

res = [ ya(1); yb(1) + 2 ];

Antes de llamar a bvp4c, debe proporcionar una conjetura para la solución que desea representar en una malla. El solver adapta después la malla mientras va ajustando la solución.

La función bvpinit construye la conjetura inicial de forma que pueda pasarse al solver bvp4c. Para una malla de [0 1 2 3 4] y conjeturas constantes de $y (x) = 1$ e $y' (x) = 0$ , la llamada a bvpinit es:

solinit = bvpinit([0 1 2 3 4],[1; 0]);

Con esta conjetura inicial, se puede resolver el problema con bvp4c. Evalúe la solución devuelta por bvp4c en algunos puntos utilizando deval y, luego, represente gráficamente los valores resultantes.

sol = bvp4c(@twoode, @twobc, solinit);

xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:));
xlabel('x')
ylabel('solution y')
hold on

Figure contains an axes object. The axes object with xlabel x, ylabel solution y contains an object of type line.

Este problema particular de valores de límites tiene exactamente dos soluciones. Puede obtener la otra solución cambiando la conjetura inicial a $y (x) = - 1$ e $y' (x) = 0$ .

solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);

xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:));
legend('Solution 1','Solution 2')
hold off

Figure contains an axes object. The axes object with xlabel x, ylabel solution y contains 2 objects of type line. These objects represent Solution 1, Solution 2.

Ecuaciones diferenciales con retardo

dde23, ddesd y ddensd resuelven ecuaciones diferenciales con retardo con varios retardos. Los ejemplos ddex1, ddex2, ddex3, ddex4 y ddex5 forman un minitutorial sobre cómo utilizar estos solvers.

El ejemplo ddex1 muestra cómo resolver el sistema de ecuaciones diferenciales

$\begin{array}{l} y_{1}^{'} (t) = y_{1} (t - 1) \\ y_{2}^{'} (t) = y_{1} (t - 1) + y_{2} (t - 0.2) \\ y_{3}^{'} (t) = y_{2} (t) . \end{array}$

Se pueden representar estas ecuaciones con la función anónima

ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];

El historial del problema (para $t \leq 0$ ) es constante:

$\begin{array}{l} y_{1} (t) = 1 \\ y_{2} (t) = 1 \\ y_{3} (t) = 1 . \end{array}$

Se puede representar el historial como un vector de unos.

ddex1hist = ones(3,1);

Un vector de dos elementos representa los retardos en el sistema de ecuaciones.

lags = [1 0.2];

Pase la función, los retardos, el historial de soluciones y el intervalo de integración $[0, 5]$ al solver como entradas. El solver genera una solución continua en todo el intervalo de integración que es adecuada para ser representada gráficamente.

sol = dde23(ddex1fun, lags, ddex1hist, [0 5]);

plot(sol.x,sol.y);
title({'An example of Wille and Baker', 'DDE with Constant Delays'});
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2','y_3','Location','NorthWest');

Figure contains an axes object. The axes object with title An example of Wille and Baker DDE with Constant Delays, xlabel time t, ylabel solution y contains 3 objects of type line. These objects represent y_1, y_2, y_3.

Ecuaciones diferenciales parciales

pdepe resuelve ecuaciones diferenciales parciales con una variable espacial y el tiempo. Los ejemplos pdex1, pdex2, pdex3, pdex4 y pdex5 forman un minitutorial sobre cómo utilizar pdepe.

Este problema de ejemplo utiliza las funciones pdex1pde, pdex1ic y pdex1bc.

pdex1pde define la ecuación diferencial

$π^{2} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} (\frac{\partial u}{\partial x}) .$

type pdex1pde

function [c,f,s] = pdex1pde(x,t,u,DuDx)
%PDEX1PDE  Evaluate the differential equations components for the PDEX1 problem.
%
%   See also PDEPE, PDEX1.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

c = pi^2;
f = DuDx;
s = 0;

pdex1ic configura la condición inicial

$u (x, 0) = \sin π x .$

type pdex1ic

function u0 = pdex1ic(x)
%PDEX1IC  Evaluate the initial conditions for the problem coded in PDEX1.
%
%   See also PDEPE, PDEX1.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

u0 = sin(pi*x);

pdex1bc configura las condiciones de límites

$u (0, t) = 0,$

$π e^{- t} + \frac{\partial}{\partial x} u (1, t) = 0 .$

type pdex1bc

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
%PDEX1BC  Evaluate the boundary conditions for the problem coded in PDEX1.
%
%   See also PDEPE, PDEX1.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

pdepe requiere la discretización espacial x y un vector de tiempos t (en el que desea obtener una instantánea de la solución). Resuelva el problema utilizando una malla de 20 nodos y solicite la solución en cinco valores de t. Extraiga y represente gráficamente el primer componente de la solución.

x = linspace(0,1,20);
t = [0 0.5 1 1.5 2];
sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

u1 = sol(:,:,1);

surf(x,t,u1);
xlabel('x');
ylabel('t');
zlabel('u');

Figure contains an axes object. The axes object with xlabel x, ylabel t contains an object of type surface.

Consulte también

ode45 | bvp4c | pdepe