此教程說明如何使用 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 和 dydtt=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')
對于較大的 μ，問題将變為剛性。此标簽表示拒絕使用普通方法計算的問題。這種情況下，要實作快速積分，需要使用特殊的數值方法。
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學習最常見的問題，每天進步一點點，業精于勤荒于嬉。
可儲存後掃碼關注哦！