Main Content

Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Ecuaciones diferenciales

Este ejemplo muestra cómo utilizar MATLAB® para formular y resolver varios tipos de ecuaciones diferenciales. MATLAB ofrece varios algoritmos numéricos para resolver una amplia gama de ecuaciones diferenciales:

  • Problemas de valor inicial

  • Problemas de valores de límites

  • Ecuaciones diferenciales con retardo

  • Ecuaciones diferenciales parciales

Problema de valor inicial

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

d2ydt2-μ(1-y2)dydt+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 de primer orden (ODE). Estas ecuaciones se evalúan para diferentes valores del parámetro μ. Para una integración más rápida, debería elegir un solver apropiado según el valor de μ.

Para μ=1, cualquiera de los solvers de ODE de MATLAB pueden resolver correctamente la ecuación de van der Pol. El solver ode45 es un ejemplo de ellos. La ecuación se resuelve en el dominio [0,20] con las condiciones iniciales y(0)=2 y dydt|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')

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

Esta solución para la ecuación de van der Pol para μ=1000 utiliza ode15s con las mismas condiciones iniciales. Debe ampliar el lapso de tiempo drásticamente 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')

Problemas de valores de límites

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

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

d2ydt2+|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ímite para el problema: y(0)=0 y 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 suposición para la solución que quiera que se represente en una red. El solver entonces adapta la red a medida que refina la solución.

La función bvpinit reúne la suposición inicial en una forma que puede pasar al solver bvp4c. Para una red de [0 1 2 3 4] y una suposición constante de y(x)=1 y y'(x)=0, la llamada a bvpinit es:

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

Con esta suposición inicial, puede resolver el problema con bvp4c. Evalúe la solución que devuelve bvp4c en algunos puntos con deval, y luego represente 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

Este problema particular de valores de límite tiene exactamente dos soluciones. Puede obtener la otra solución si cambia la suposición inicial a y(x)=-1 y 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

Ecuaciones diferenciales con retardo

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

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

y1'(t)=y1(t-1)y2'(t)=y1(t-1)+y2(t-0.2)y3'(t)=y2(t).

Puede 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 t0) es constante:

y1(t)=1y2(t)=1y3(t)=1.

Puede representar el historial como un vector de unos.

ddex1hist = ones(3,1);

Un vector de dos elementos representa los retardos del 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 produce una solución continua sobre todo el intervalo de integración que sea susceptible de representar.

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');

Ecuaciones diferenciales parciales

pdepe resuelve ecuaciones diferenciales parciales en una variable de espacio y tiempo. Los ejemplos pdex1, pdex2, pdex3, pdex4 y pdex5 componen un minitutorial sobre cómo utilizar pdepe.

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

pdex1pde define la ecuación diferencial

π2ut=x(ux).

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 establece 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 establece las condiciones de límites

u(0,t)=0,

πe-t+xu(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 exige una discretización espacial x y un vector de tiempo t (en el que le interesa una instantánea de la solución). Resuelva el problema con una red de 20 nodos y solicite la solución en cinco valores de t. Extraiga y represente 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');

Consulte también

| |

Temas relacionados