此教程说明如何使用 MATLAB 构造几种不同类型的微分方程并求解。MATLAB 提供了多种数值算法来求解各种微分方程：

• 初始值问题

• 边界值问题

• 时滞微分方程

• 偏微分方程

初始值问题

vanderpoldemo 是用于定义 van der Pol 方程的函数

type vanderpoldemo

function dydt = vanderpoldemo(t,y,Mu)dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

该方程写作包含两个一阶常微分方程 (ODE) 的方程组。将针对参数 μ 的不同值计算这些方程。为了实现更快的积分，您应该根据 μ 的值选择合适的求解器。

当 μ=1 时，任何 MATLAB ODE 求解器都能有效地求解 van der Pol 方程。ode45 求解器就是其中之一。该方程在域 [0,20] 中求解，初始条件为 y(0)=2 和 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 solutionplot(t,y(:,1))xlabel('t')ylabel('solution y')title('van der Pol Equation, \mu = 1')

对于较大的 μ，问题将变为刚性。此标签表示拒绝使用普通方法计算的问题。这种情况下，要实现快速积分，需要使用特殊的数值方法。ode15s、ode23s、ode23t 和 ode23tb 函数可有效地求解刚性问题。

当 μ=1000 时，van der Pol 方程的求解使用 ode15s，初始条件相同。您需要将时间范围大幅度延长到 [0,3000] 才能看到解的周期性变化。

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.dydx = [ y(2); -abs(y(1)) ];

函数 twobc 求解该问题的边界条件为：y(0)=0 和 y(4)=−2。

type twobc

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

在调用 bvp4c 之前，您必须为要在网格中表示的解提供一个猜想值。然后，求解器就像对解进行平滑处理一样修改网格。

bvpinit 函数以您可以传递给求解器 bvp4c 的形式设定初始猜想值。对于 [0 1 2 3 4] 的网格以及 y(x)=1 和 y'(x)=0 的常量猜想值，对 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

此特定的边界值问题实际上有两种解。通过将初始猜想值更改为 y(x)=−1 和 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

时滞微分方程

     dde23、ddesd 和 ddensd 可以求解具有各种时滞的时滞微分方程。示例 ddex1、ddex2、ddex3、ddex4 和 ddex5 构成了这些求解器的迷你使用教程。

ddex1 示例说明如何求解微分方程组

y′1(t)=y1(t−1)y′2(t)=y1(t−1)+y2(t−0.2)y′3(t)=y2(t).

您可以使用匿名函数表示这些方程

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

问题的历史解(t≤0 时)固定不变：

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

您可以将历史解表示为由 1 组成的向量。

ddex1hist = ones(3,1);

采用二元素向量表示方程组中的时滞。

lags = [1 0.2];

将函数、时滞、历史解和积分区间 [0,5] 作为输入传递给求解器。求解器在整个积分区间生成适合绘图的连续解。

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.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.u0 = sin(pi*x);

pdex1bc 设置边界条件

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

关注公众号：MATLAB基于模型的设计 (ID：xaxymaker) ，每天推送MATLAB学习最常见的问题，每天进步一点点，业精于勤荒于嬉。

可保存后扫码关注哦！