目錄

理論知識(shí)

一、概念

二、解法

matlab微分方程求解

一、解析解

1.1 解析解的存在

1.2 解析解的解法

1.3 實(shí)例

二、數(shù)值解

2.1 概述

2.2 優(yōu)化措施

2.3 解法

2.4?檢驗(yàn)

---

理論知識(shí)

一、概念

微分方程：含導(dǎo)數(shù)或微分的方程。

解：滿足微分方程的函數(shù)。

特解/通解：特解指的是滿足微分方程的某一個(gè)解；通解指的是滿足微分方程的一組解。

階：微分方程中導(dǎo)數(shù)或微分的最高階數(shù)。

線性/非線性：（幾何意義：疊加原理）方程中的函數(shù)和它的各階導(dǎo)數(shù)都是一次方為線性微分方程，否則為非線性。例：

• y'=sin(x)*y 線性
• y'=y^2 非線性

齊次/非齊次：（代數(shù)意義：次數(shù)）齊次微分方程中不含常數(shù)項(xiàng),也不含僅由x的各種運(yùn)算組合構(gòu)成的項(xiàng)（比如4xx,sinx等），否則為非齊次。

常微分/偏微分：未知函數(shù)是一元函數(shù)的，叫常微分方程；未知函數(shù)是多元函數(shù)的叫做偏微分方程。

初值問題和邊值問題：附加條件中未知函數(shù)及其導(dǎo)數(shù)的獨(dú)立變量取值相同，則為初值問題；附加條件中未知函數(shù)及其導(dǎo)數(shù)的獨(dú)立變量取值不同，則為邊界值問題。例：

• y(0)=1,y'(0)=2 初值問題
• y(0)=1,y'(1)=2?邊值問題

二、解法

在matlab中解微分方程，方法有兩種：

解析解方法：嚴(yán)格按照公式邏輯推導(dǎo)得到的，具有基本的函數(shù)形式。

數(shù)值解方法：采用某種計(jì)算方法，在特定的條件下得到的一個(gè)近似數(shù)值結(jié)果，如有限元法，數(shù)值逼近法，插值法等等。

解析解方法可回顧《高等數(shù)學(xué)》中相關(guān)公式。

數(shù)值解方法可回顧《數(shù)值分析》中相關(guān)解法。

matlab微分方程求解

一、解析解

1.1 解析解的存在

由Abel-Ruffini定理，四次及以下的多項(xiàng)式代數(shù)方程是能求出根的解析解的，即低階常系數(shù)線性微分方程有一般意義下的解析解。

非線性微分方程只能用數(shù)值解法求解，即使看起來很簡(jiǎn)單的非線性微分方程也是沒有解析解的，只有極特殊的非線性微分方程解析可解。

1.2 解析解的解法

利用dsolve函數(shù)

• S = dsolve(eqn)
• S = dsolve(eqn,cond)
• S = dsolve(___,Name,Value)
• [y1,...,yN] = dsolve(___)

1.3 實(shí)例

例1.1?輸入信號(hào)為u(t)=exp(-5*t)*cos(2*t+1)+5，求微分方程diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y=5*diff(u,t,2)+4*diff(u,t)+2*u的通解。初值條件：y(0)=3,y1(0)=2,y2(0)=0,y3(0)=0，求方程的特解。（約定y1代表一階導(dǎo)數(shù)，以此類推）

代碼如下：

% y=dsolve(f1,f2,...,fm) 默認(rèn)自變量為t % y=dsolve(f1,f2,...,fm,'x') 指明自變量 % f可由字符串表示也可由符號(hào)表達(dá)式表示%用字符串表達(dá)式，新版本會(huì)被移除 syms t; u=exp(-5*t)*cos(2*t+1)+5; uu=5*diff(u,t,2)+4*diff(u,t)+2*u; y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',char(uu)]); %求解 y=simplify(y) %化簡(jiǎn)%用符號(hào)表達(dá)式，推薦 syms y(t); eqn=diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y==uu; y=dsolve(eqn); y=simplify(y) %檢驗(yàn) diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y-uu%求特解并繪圖 syms y(t); eqn=diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y==uu; y1=diff(y);y2=diff(y,2);y3=diff(y,3);y4=diff(y,4); % 需要引入中間變量 % 用字符串求解的情況，不需要引入中間變量。'y(0)=3','Dy(0)=2','D2y(0)=0'... cond=[y(0)==3,y1(0)==2,y2(0)==0,y3(0)==0]; z=dsolve(eqn,cond); z=simplify(z) ezplot(z,[0,5]) fplot(z,[0,5]) % ezplot(fun,[xmin,xmax])繪制fun(x)在以下域上的圖形:xmin<x<xmax % ezplot適用隱式函數(shù) % 推薦fplot代替ezplot double(subs(z,t,5)) % 驗(yàn)證圖像，對(duì)自變量賦值，并求小數(shù)解% 改變初值條件，再求特解，用字符串表示方式 z1=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',char(uu)],...'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5') %求解 fplot(z1,[0,5])

運(yùn)行結(jié)果：

y =
C1*exp(-4*t) - (547*exp(-5*t)*sin(2*t + 1))/520 - (343*exp(-5*t)*cos(2*t + 1))/520 + C2*exp(-3*t) + C3*exp(-2*t) + C4*exp(-t) + 5/12

ans =0

z =
19*exp(-t) - (69*exp(-2*t))/2 + (73*exp(-3*t))/3 - (25*exp(-4*t))/4 + (97*exp(-t)*sin(1))/60 - (51*exp(-2*t)*sin(1))/13 + (5*exp(-3*t)*sin(1))/8......

ans =?0.5679?

z1 =?
(exp(-5*t)*(445*cos(2*t + 1) - 65*exp(5*t) + 102*sin(2*t + 1)))/26 - (exp(-5*t)*(537*cos(2*t + 1) - 40*exp(5*t) + 15*sin(2*t + 1)))/24 - (exp(-5*t)*(25*exp(5*t) - 542*cos(2*t + 1) ......

例1.2?求解有復(fù)數(shù)極點(diǎn)的微分方程diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y=diff(u)+3*u，輸入信號(hào)：u(t)=sin(t)，初值條件：y(0)=0,y1(0)=0,y2(0)=0,y3(0)=0,y4(0)=0

代碼如下：

% 有復(fù)數(shù)極點(diǎn)的微分方程 syms y(t) u(t); u=sin(t); uu=diff(u)+3*u; eqn=diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y==uu; y1=diff(y);y2=diff(y,2);y3=diff(y,3);y4=diff(y,4); % 需要引入中間變量 cond=[y(0)==0,y1(0)==0,y2(0)==0,y3(0)==0,y4(0)==0]; y=dsolve(eqn,cond); y=simplify(y) fplot(y,[0,30]) % 檢驗(yàn) diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y-uu

?運(yùn)行結(jié)果：

y =
exp(-t) - cos(t)/5 - (2*sin(t))/5 - (4*exp(-t)*cos(t))/5 + (11*exp(-t)*sin(t))/10 - (t*exp(-t)*cos(t))/2

ans =0

?例1.3?求線性微分方程組

• diff(y)==4*x+3*y+4*exp(-t))
• diff(x,2)+2*diff(x)==x+2*y-exp(-t)

代碼如下：

%微分方程組,無法檢驗(yàn)結(jié)果 syms y(t) x(t); eqn1=diff(x,2)+2*diff(x)==x+2*y-exp(-t); eqn2=diff(y)==4*x+3*y+4*exp(-t); [x,y]=dsolve(eqn1,eqn2) diff(x,2)+2*diff(x)-x+2*y-exp(-t) diff(y)-4*x+3*y+4*exp(-t)

?運(yùn)行結(jié)果：

x =
exp(t*(6^(1/2) + 1))*(6^(1/2)/5 - 1/5)*(C2 + exp(- 2*t - 6^(1/2)*t)*((11*6^(1/2))/3 - 37/4)) - exp(-t)*(C1 + 6*t) - exp(-t*(6^(1/2) - 1))*(6^(1/2)/5 + 1/5)*(C3 - exp(6^(1/2)*t - 2*t)*((11*6^(1/2))/3 + 37/4))
?
y =
exp(-t)*(C1 + 6*t) + exp(t*(6^(1/2) + 1))*((2*6^(1/2))/5 + 8/5)*(C2 + exp(- 2*t - 6^(1/2)*t)*((11*6^(1/2))/3 - 37/4)) - exp(-t*(6^(1/2) - 1))*((2*6^(1/2))/5 - 8/5)*(C3 - exp(6^(1/2)*t - 2*t)*((11*6^(1/2))/3 + 37/4))
?
ans =(無法檢驗(yàn))
4*exp(-t)*(C1 + 6*t) - exp(-t) - exp(t*(6^(1/2) + 1))*(6^(1/2)/5 - 1/5)*(C2 + exp(- 2*t - 6^(1/2)*t)*((11*6^(1/2))/3 - 37/4)) + ......
?
ans =(無法檢驗(yàn))
10*exp(-t) + 6*exp(-t)*(C1 + 6*t) - 4*exp(t*(6^(1/2) + 1))*(6^(1/2)/5 - 1/5)*(C2 + exp(- 2*t - 6^(1/2)*t)*((11*6^(1/2))/3 - 37/4)) +......

?例1.4?求解自變量為x的時(shí)變線性微分方程(2*x+3)^3*diff(y,3)+3*(2*x+3)*diff(y)-6*y=0

代碼如下：

% 自變量為x的時(shí)變微分方程 syms x y(x); y=dsolve((2*x+3)^3*diff(y,3)+3*(2*x+3)*diff(y)-6*y==0)； simplify(y) (2*x+3)^3*diff(y,3)+3*(2*x+3)*diff(y)-6*y

運(yùn)行結(jié)果：

y =
C1*(x + 3/2) - 2*C2*(x + 3/2)^(1/2) + C3*(2*x + 6)*(x + 3/2)^(1/2)

ans =(無法檢驗(yàn))
12*C2*(x + 3/2)^(1/2) - 6*C1*(x + 3/2) - (2*x + 3)^3*((3*C3)/(2*(x + 3/2)^(3/2)) + (3*C2)/(4*(x + 3/2)^(5/2)) - (3*C3*(2*x + 6))/(8*(x + 3/2)^(5/2))) + (6*x + 9)*(C1 - C2/(x + 3/2)^(1/2) + 2*C3*(x + 3/2)^(1/2) + (C3*(2*x + 6))/(2*(x + 3/2)^(1/2))) - 6*C3*(2*x + 6)*(x + 3/2)^(1/2)
?

例1.5?求解高階常系數(shù)微分方程組

• diff(x,2)-x+y+z=0
• x+diff(y,2)-y+z=0
• +y+diff(z,2)-z=0

?代碼如下：

% 高階常系數(shù)線性微分方程 syms x(t) y(t) z(t); eqn1=diff(x,2)-x+y+z==0; eqn2=x+diff(y,2)-y+z==0; eqn3=x+y+diff(z,2)-z==0; x1=diff(x);y1=diff(y);z1=diff(z); cond=[x(0)==1,y(0)==0,z(0)==0,x1(0)==0,y1(0)==0,z1(0)==0]; [x y z]=dsolve(eqn1,eqn2,eqn3,cond) % 檢驗(yàn)，成功 diff(x,2)-x+y+z x+diff(y,2)-y+z x+y+diff(z,2)-z

運(yùn)行結(jié)果：

x =
exp(2^(1/2)*t)/3 + exp(-2^(1/2)*t)/3 + cos(t)/3
?
y =
cos(t)/3 - exp(-2^(1/2)*t)/6 - exp(2^(1/2)*t)/6
?
z =
cos(t)/3 - exp(-2^(1/2)*t)/6 - exp(2^(1/2)*t)/6
?
ans =0?

ans =0?

ans =0

例1.6?求解時(shí)變線性微分方程x^2*(2*x-1)*diff(y,x,3)+(4*x-3)*x*diff(y,x,2)-2*x*diff(y,x)+2*y=0

代碼如下：

% 自變量為x的時(shí)變微分方程 syms x y(x); eqn=x^2*(2*x-1)*diff(y,x,3)+(4*x-3)*x*diff(y,x,2)-2*x*diff(y,x)+2*y==0; y=dsolve(eqn)； simplify(y) % 檢驗(yàn)，失敗 x^2*(2*x-1)*diff(y,3)-(4*x-3)*x*diff(y,2)-2*x*diff(y)+2*y

運(yùn)行結(jié)果：

y =
-(2*C1 - C3 + 8*C3*x + 32*C2*x^2 + 8*C3*x^2*log(x))/(16*x)

ans =(無法檢驗(yàn))
2*x*((8*C3 + 64*C2*x + 8*C3*x + 16*C3*x*log(x))/(16*x) - (2*C1 - C3 + 8*C3*x + 32*C2*x^2 + 8*C3*x^2*log(x))/(16*x^2))+...
?

例1.7?求解特殊非線性微分方程diff(x)=x*(1-x^2)

代碼如下：

% 特殊非線性方程 syms x(t); x=dsolve(diff(x)==x*(1-x^2))

運(yùn)行結(jié)果：

x =
?(-1/(exp(C1 - 2*t) - 1))^(1/2)
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1
? ? ? ? ? ? ? ? ? ? ? ? ? ? ?-1

注：若改為x=dsolve(diff(x)==x*(1-x^2)+1)，則無解。

二、數(shù)值解

2.1 概述

大部分為關(guān)于微分方程的初值問題的數(shù)值解法，這類問題需要用一階顯式微分方程組描述為

x ' ( t ) = f ( t , x ( t ) )

• x(t)=[x1(t);x2(t);...xn(t)] 狀態(tài)向量
• f(.)=[f1(.);f2(.);...fn(.)]?任意非線性函數(shù)
• x0=[x1(t0);x2(t0);...xn(t0)]?初始狀態(tài)
• [t0,tn]?時(shí)間區(qū)間（tn為終止時(shí)間）?

2.2 優(yōu)化措施

• 選擇適當(dāng)?shù)牟介L(zhǎng)h：不能太大和太小，太大精度不夠，太小減慢計(jì)算精度，增加計(jì)算次數(shù)從而增加累積誤差。
• 改進(jìn)近似算法精度：Euler算法將原始積分進(jìn)行梯形近似，精度低；精度較高的是Runge-Kutta法和Adams法。
• 采用變步長(zhǎng)方法：采用自適應(yīng)步長(zhǎng)算法。

2.3 解法

1、手動(dòng)編寫算法。例如：

四階定步長(zhǎng)Runge-Kutta算法matlab代碼如下：

function [tout,yout]=rk_4(fun,tspan,y0) ts=tspan; t0=ts(1); tf=ts(2); yout=[]; tout=[]; y0=y0(:); if length(tspan)==3,h=ts(3); else, h=(ts(2)-ts(1))/100;tf=ts(2); end for t = [t0:h:tf] % 對(duì)各個(gè)時(shí)間點(diǎn)循環(huán)計(jì)算k1=h*fun(t,y0);k2=h*fun(t+h/2,y0+k1/2);k3=h*fun(t+h/2,y0+k2/2);k4=h*fun(t+h,y0+k3);y0=y0+(k1+2*k2+2*k3+k4)/6;yout=[yout;y0.'];tout=[tout;t]; end

2、ode系列函數(shù)：最常用ode45()

介紹

ode45()實(shí)現(xiàn)了變步長(zhǎng)四階五級(jí)?Runge-Kutta-Felhberg算法，可以使用變步長(zhǎng)的方法求解微分方程。ode系列函數(shù)功能總結(jié)如下：

用法

函數(shù)調(diào)用：

• [t,y] = ode45(odefun,tspan,y0)
• [t,y] = ode45(odefun,tspan,y0,options)
• [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)
• sol = ode45(___)

微分方程/組描述：

• function xd=funname(t,x) % 不需要附加參數(shù)格式
• function xd=funname(t,x,p1,p2,...,pm) % 可以使用附加參數(shù)

控制條件：使用odesrt()函數(shù)創(chuàng)建或修改 options 結(jié)構(gòu)體。

• options=odeset('RelTol',1e-7)
• options=odeset;options.RelTol=1e-7

例：options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) 指定 1e-5 的相對(duì)誤差容限、打開求解器統(tǒng)計(jì)信息的顯示，并指定輸出函數(shù) @odeplot 在計(jì)算時(shí)繪制解。

實(shí)例

例2.1?解微分方程組

• diff(x1)=-beta*x1+x2*x3
• diff(x2)=-Pho*x2+Pho*x3
• diff(x3)=-x1*x2+sigma*x2-x3

?其中，beta=8/3,Pho=10,sigma=28，初值:x1(0)=x2(0)=0,x3(0)=10e-10

代碼1：（不帶有附加參數(shù)的函數(shù)）

f=@(t,x) [-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)] % 編寫匿名函數(shù)描述動(dòng)態(tài)模型 tn=100; x0=[0;0;1e-10]; [t,x]=ode45(f,[0,tn],x0); %求方程數(shù)值解 plot(t,x); figure; plot3(x(:,1),x(:,2),x(:,3)); %相空間軌跡三維圖 grid % 顯示或隱藏坐標(biāo)區(qū)網(wǎng)格線 figure; comet3(x(:,1),x(:,2),x(:,3)); %繪制動(dòng)畫式軌跡

運(yùn)行結(jié)果：

代碼2：（帶有附加參數(shù)的函數(shù)）

函數(shù)程序：

function f=c_7_2_fun % 定義一個(gè)接口，通過匿名函數(shù)調(diào)用 % 通過定義接口的方式，在一個(gè)函數(shù)文件實(shí)例里可編寫多個(gè)函數(shù) f.lorenz1=@lorenz1; endfunction dx=lorenz1(t,x,b,r,s) % 帶有附加參數(shù)的微分方程描述曲線 dx=[-b*x(1)+x(2)*x(3);-r*x(2)+r*x(3);-x(1)*x(2)+s*x(2)-x(3)]; end

主程序：

%% 編寫帶有附加參數(shù)的函數(shù) % 函數(shù)需列在單獨(dú)的文件里，且函數(shù)名和文件名一致，可通過定義接口的方式，在一個(gè)函數(shù)文件例里編寫多個(gè)函數(shù) f=c_7_2_fun % 調(diào)用接口 % 設(shè)置附加參數(shù)的值 b1=8/3;r1=10;s1=28;tn=100;x0=[0;0;1e-10]; [t,x]=ode45(@f.lorenz1,[0,tn],x0,[],b1,r1,s1); % 無需使用和函數(shù)一樣的變量名 % 求方程數(shù)值解 % 空矩陣表示使用默認(rèn)的控制模版 figure(1); plot(t,x); figure(2); comet3(x(:,1),x(:,2),x(:,3)); %繪制動(dòng)畫式軌跡 grid%% 改變附加參數(shù)的值 b2=2;r2=5;s2=20;tn=100;x0=[0;0;1e-10]; [t2,x2]=ode45(@f.lorenz1,[0,tn],x0,[],b2,r2,s2); % 求方程數(shù)值解 % 空矩陣表示使用默認(rèn)的控制模版 figure(3); plot(t2,x2); figure(4); comet3(x2(:,1),x2(:,2),x2(:,3)); %繪制動(dòng)畫式軌跡 grid %只針對(duì)上一個(gè)繪圖% TIPS %% 使用匿名函數(shù)則可以不用附加參數(shù) b=8/3;r=10;s=28;tn=100;x0=[0;0;1e-10]; f=@(t,x) [-b*x(1)+x(2)*x(3);-r*x(2)+r*x(3);-x(1)*x(2)+s*x(2)-x(3)]; % 匿名函數(shù)可以直接使用matlab工作空間中的變量 [t,x]=ode45(f,[0,tn],x0); %求方程數(shù)值解 plot(t,x); figure; plot3(x(:,1),x(:,2),x(:,3)); %相空間軌跡三維圖 grid % 顯示或隱藏坐標(biāo)區(qū)網(wǎng)格線

注：函數(shù)需列在單獨(dú)的文件里，且函數(shù)名和文件名一致，才可保證正常調(diào)用，即一個(gè)函數(shù)對(duì)應(yīng)一個(gè)函數(shù)文件；現(xiàn)可通過定義接口的方式，在一個(gè)函數(shù)文件實(shí)例里編寫多個(gè)函數(shù)，調(diào)用時(shí)調(diào)用接口。

運(yùn)行結(jié)果：（改變參數(shù)前的運(yùn)行結(jié)果同上）

2.4?檢驗(yàn)

• 修改仿真控制參數(shù)，例如：將‘RelTol’或'AbsTol'選項(xiàng)設(shè)置成一個(gè)更小的值，觀察結(jié)果是否和上次一致，如有不可接受的差異，應(yīng)考慮進(jìn)一步減小誤差值。
• 選擇不同的微分方程求解算法。

本文僅涉及一般常微分方程求解，對(duì)于一些特殊常微分方程的求解方法見同專欄相關(guān)文章。

總結(jié)

以上是生活随笔為你收集整理的matlab求解常微分方程——从原理到实践（代码详解）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。