Ecuaciones diferenciales
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 amplia variedad 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
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 (EDO) 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 , cualquiera de los solvers de EDO 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 con las condiciones iniciales y .
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 convierte en rígido. Esta etiqueta se aplica a problemas que resisten los intentos a 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 para la ecuación de van der Pol para utiliza ode15s
con las mismas condiciones iniciales. Se necesita ampliar drásticamente el periodo de tiempo a 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 para ecuaciones diferenciales ordinarias.
La función de ejemplo twoode
tiene una ecuación diferencial escrita como un sistema de dos EDO de primer orden. La ecuación diferencial es
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 .
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 , 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
Este problema particular de valores de límites tiene exactamente dos soluciones. Puede obtener la otra solución cambiando la conjetura inicial a y .
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
forman un minitutorial sobre como utilizar estos solvers.
El ejemplo ddex1
muestra cómo resolver el sistema de ecuaciones diferenciales
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 ) es constante:
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 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');
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
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
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
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');