微分方程
此示例说明如何使用 MATLAB® 构造几种不同类型的微分方程并求解。MATLAB 提供了多种数值算法来求解各种微分方程:
初始值问题
边界值问题
时滞微分方程
偏微分方程
初始值问题
vanderpoldemo
是用于定义 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)];
该方程写作包含两个一阶常微分方程 (ODE) 的方程组。将针对参数 的不同值计算这些方程。为了实现更快的积分,您应该根据 的值选择合适的求解器。
当 时,任何 MATLAB ODE 求解器都能有效地求解 van der Pol 方程。ode45
求解器就是其中之一。该方程在域 中求解,初始条件为 和 。
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')
对于较大的 ,问题将变为刚性。此标签表示拒绝使用普通方法计算的问题。这种情况下,要实现快速积分,需要使用特殊的数值方法。ode15s
、ode23s
、ode23t
和 ode23tb
函数可有效地求解刚性问题。
当 时,van der Pol 方程的求解使用 ode15s
,初始条件相同。您需要将时间范围大幅度延长到 才能看到解的周期性变化。
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')
边界值问题
bvp4c
和 bvp5c
可以求解常微分方程的边界值问题。
示例函数 twoode
将一个微分方程写作包含两个一阶 ODE 的方程组。此微分方程为
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)) ];
函数 twobc
求解该问题的边界条件为: 和 。
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 ];
在调用 bvp4c
之前,您必须为要在网格中表示的解提供一个猜想值。然后,求解器就像对解进行平滑处理一样修改网格。
bvpinit
函数以您可以传递给求解器 bvp4c
的形式设定初始猜想值。对于 [0 1 2 3 4]
的网格以及 和 的常量猜想值,对 bvpinit
的调用为:
solinit = bvpinit([0 1 2 3 4],[1; 0]);
利用这个初始猜想值,您可以使用 bvp4c
对该问题求解。使用 deval
计算 bvp4c
在某些点返回的解,然后绘制结果值。
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
此特定的边界值问题实际上有两种解。通过将初始猜想值更改为 和 ,可以求出另一个解。
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
时滞微分方程
dde23
、ddesd
和 ddensd
可以求解具有各种时滞的时滞微分方程。示例 ddex1
、ddex2
、ddex3
、ddex4
和 ddex5
构成了这些求解器的迷你使用教程。
ddex1
示例说明如何求解微分方程组
您可以使用匿名函数表示这些方程
ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];
问题的历史解( 时)固定不变:
您可以将历史解表示为由 1 组成的向量。
ddex1hist = ones(3,1);
采用二元素向量表示方程组中的时滞。
lags = [1 0.2];
将函数、时滞、历史解和积分区间 作为输入传递给求解器。求解器在整个积分区间生成适合绘图的连续解。
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');
偏微分方程
pdepe
使用一个空间变量和时间对偏微分方程求解。示例 pdex1
、pdex2
、pdex3
、pdex4
和 pdex5
构成了 pdepe
的迷你使用教程。
此示例问题使用函数 pdex1pde
、pdex1ic
和 pdex1bc
。
pdex1pde
定义微分方程
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
设置初始条件
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
设置边界条件
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
需要提供空间离散 x
和时间向量 t
(您要获取解快照的时间点)。使用包含 20 个节点的网格求解此问题,并请求五个 t
值的解。提取解的第一个分量并绘图。
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');