目錄

理論知識

一、概念

二、解法

matlab微分方程求解

一、解析解

1.1 解析解的存在

1.2 解析解的解法

1.3 實例

二、數值解

2.1 概述

2.2 優化措施

2.3 解法

2.4?檢驗

---

理論知識

一、概念

微分方程：含導數或微分的方程。

解：滿足微分方程的函數。

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

階：微分方程中導數或微分的最高階數。

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

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

齊次/非齊次：（代數意義：次數）齊次微分方程中不含常數項,也不含僅由x的各種運算組合構成的項（比如4xx,sinx等），否則為非齊次。

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

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

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

二、解法

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

解析解方法：嚴格按照公式邏輯推導得到的，具有基本的函數形式。

數值解方法：采用某種計算方法，在特定的條件下得到的一個近似數值結果，如有限元法，數值逼近法，插值法等等。

解析解方法可回顧《高等數學》中相關公式。

數值解方法可回顧《數值分析》中相關解法。

matlab微分方程求解

一、解析解

1.1 解析解的存在

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

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

1.2 解析解的解法

利用dsolve函數

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

1.3 實例

例1.1?輸入信號為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代表一階導數，以此類推）

代碼如下：

% y=dsolve(f1,f2,...,fm) 默認自變量為t % y=dsolve(f1,f2,...,fm,'x') 指明自變量 % f可由字符串表示也可由符號表達式表示%用字符串表達式，新版本會被移除 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) %化簡%用符號表達式，推薦 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) %檢驗 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適用隱式函數 % 推薦fplot代替ezplot double(subs(z,t,5)) % 驗證圖像，對自變量賦值，并求小數解% 改變初值條件，再求特解，用字符串表示方式 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 =
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?求解有復數極點的微分方程diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y=diff(u)+3*u，輸入信號：u(t)=sin(t)，初值條件：y(0)=0,y1(0)=0,y2(0)=0,y3(0)=0,y4(0)=0

代碼如下：

% 有復數極點的微分方程 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]) % 檢驗 diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y-uu

?運行結果：

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)

代碼如下：

%微分方程組,無法檢驗結果 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)

?運行結果：

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 =(無法檢驗)
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 =(無法檢驗)
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的時變線性微分方程(2*x+3)^3*diff(y,3)+3*(2*x+3)*diff(y)-6*y=0

代碼如下：

% 自變量為x的時變微分方程 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 =
C1*(x + 3/2) - 2*C2*(x + 3/2)^(1/2) + C3*(2*x + 6)*(x + 3/2)^(1/2)

ans =(無法檢驗)
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?求解高階常系數微分方程組

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

?代碼如下：

% 高階常系數線性微分方程 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) % 檢驗，成功 diff(x,2)-x+y+z x+diff(y,2)-y+z x+y+diff(z,2)-z

運行結果：

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?求解時變線性微分方程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的時變微分方程 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) % 檢驗，失敗 x^2*(2*x-1)*diff(y,3)-(4*x-3)*x*diff(y,2)-2*x*diff(y)+2*y

運行結果：

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

ans =(無法檢驗)
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))

運行結果：

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

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

二、數值解

2.1 概述

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

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

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

2.2 優化措施

• 選擇適當的步長h：不能太大和太小，太大精度不夠，太小減慢計算精度，增加計算次數從而增加累積誤差。
• 改進近似算法精度：Euler算法將原始積分進行梯形近似，精度低；精度較高的是Runge-Kutta法和Adams法。
• 采用變步長方法：采用自適應步長算法。

2.3 解法

1、手動編寫算法。例如：

四階定步長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] % 對各個時間點循環計算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系列函數：最常用ode45()

介紹

ode45()實現了變步長四階五級?Runge-Kutta-Felhberg算法，可以使用變步長的方法求解微分方程。ode系列函數功能總結如下：

用法

函數調用：

• [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) % 不需要附加參數格式
• function xd=funname(t,x,p1,p2,...,pm) % 可以使用附加參數

控制條件：使用odesrt()函數創建或修改 options 結構體。

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

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

實例

例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：（不帶有附加參數的函數）

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)] % 編寫匿名函數描述動態模型 tn=100; x0=[0;0;1e-10]; [t,x]=ode45(f,[0,tn],x0); %求方程數值解 plot(t,x); figure; plot3(x(:,1),x(:,2),x(:,3)); %相空間軌跡三維圖 grid % 顯示或隱藏坐標區網格線 figure; comet3(x(:,1),x(:,2),x(:,3)); %繪制動畫式軌跡

運行結果：

代碼2：（帶有附加參數的函數）

函數程序：

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

主程序：

%% 編寫帶有附加參數的函數 % 函數需列在單獨的文件里，且函數名和文件名一致，可通過定義接口的方式，在一個函數文件例里編寫多個函數 f=c_7_2_fun % 調用接口 % 設置附加參數的值 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); % 無需使用和函數一樣的變量名 % 求方程數值解 % 空矩陣表示使用默認的控制模版 figure(1); plot(t,x); figure(2); comet3(x(:,1),x(:,2),x(:,3)); %繪制動畫式軌跡 grid%% 改變附加參數的值 b2=2;r2=5;s2=20;tn=100;x0=[0;0;1e-10]; [t2,x2]=ode45(@f.lorenz1,[0,tn],x0,[],b2,r2,s2); % 求方程數值解 % 空矩陣表示使用默認的控制模版 figure(3); plot(t2,x2); figure(4); comet3(x2(:,1),x2(:,2),x2(:,3)); %繪制動畫式軌跡 grid %只針對上一個繪圖% TIPS %% 使用匿名函數則可以不用附加參數 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)]; % 匿名函數可以直接使用matlab工作空間中的變量 [t,x]=ode45(f,[0,tn],x0); %求方程數值解 plot(t,x); figure; plot3(x(:,1),x(:,2),x(:,3)); %相空間軌跡三維圖 grid % 顯示或隱藏坐標區網格線

注：函數需列在單獨的文件里，且函數名和文件名一致，才可保證正常調用，即一個函數對應一個函數文件；現可通過定義接口的方式，在一個函數文件實例里編寫多個函數，調用時調用接口。

運行結果：（改變參數前的運行結果同上）

2.4?檢驗

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

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。