目錄

微分方程的轉換

一、單個高階常微分方程

二、高階常微分方程組

剛性微分方程求解

隱式微分方程求解

微分代數方程求解

---

微分方程的轉換

根據微分方程求解的標準型，要得到微分方程的數值解，應該先將該方程變換成一階常微分方程組。

一、單個高階常微分方程

1、原狀態

高階常微分方程一般形式：

初值為y(0),y'(0),...,y^(n-1)(0)

2、變換后

可以選擇一組狀態變量x1=y,x2=y',...,xn=y^(n-1)，可以將原方程變換為一階方程組形式：

• x'1=x2
• x'2=x3
• ...
• x'n=f(t,x1,x2,...,xn)

初值為x1(0)=y(0),x2(0)=y'(0),...,xn(0)=y^(n-1)(0)

例1.1?求y''+mu*(y^2-1)*y'+y=0，初值為y(0)=0.2,y'(0)=-0.7

代碼如下：

%% 單個高階常微分方程處理方法 f=@(t,x,mu) [x(2);-mu*(x(1)^2-1)*x(2)-x(1)];% 帶附加參數的匿名函數 x0=[-0.2;-0.7]; tn=20;mu=1; tic; % tic，toc記錄時間 [t1,y1]=ode45(f,[0,tn],x0,[],mu); toc; mu=2; tic; [t2,y2]=ode45(f,[0,tn],x0,[],mu); toc; figure(1); plot(t1,y1(:,1),t2,y2(:,1),'--'); figure(2); plot(t1,y1(:,2),t2,y2(:,2),'--'); figure(3); plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),'--'); % 繪制相平面圖 % mu=1000，tn=3000，x0=[2;0]時,未出現內存不夠的情況,但用求解剛性微分方程的方式計算量小得多 % options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);是解微分方程時的選項設置 % 'RelTol',1e-4,是相對誤差設置，'AbsTol',[1e-4 1e-4 1e-5]是絕對誤差設置 f=@(t,x,mu) [x(2);-mu*(x(1)^2-1)*x(2)-x(1)]; x0=[2;0];tn=3000; mu=1000; tic; options=odeset('RelTol',1e-10,'AbsTol',1e-10); [t3,y3]=ode45(f,[0,tn],x0,[],mu);toc;figure(4); plot(t3,y3(:,1)); figure(5); plot(t3,y3(:,2),'--'); figure(6); plot(y3(:,1),y3(:,2)); %繪制相平面圖

運行結果：

mu=1和mu=2時，實線為[t,y]；虛線為[t,y']；用時分別為?0.002099 秒和0.001892 秒。

??mu=1000，tn=3000，x0=[2;0]時，實線為[t,y]；虛線為[t,y']；用時為28.018461 秒。

注：用ode45()計算mu=1000，tn=3000，x0=[2;0]時的方程，如果電腦內存不夠的話會報錯，因為步長過小導致，可采取剛性方程的算法，用ode15s可大大減小計算量。

二、高階常微分方程組

1、原狀態

如果可以顯式地將兩個方程寫成：

初值為y(0),y'(0),...,y^(n-1)(0)

2、變換后

仍然可以選擇一組狀態變量x1=x,x2=x',...,xm=x^(m-1)，x(m+1)=y,x(m+2)=y',...,x(m+n)=y^(n-1)，可以將原方程變換為一階微分方程組：

• x'1=x2
• ...
• x'm=f(t,x1,x2,...,x(m+n))
• x'(m+1)=x(m+2)
• ...
• x'(m+n)=g(t,x1,x2,...,x(m+n))

例1.2?衛星的運動軌跡滿足下面方程組：

• x''=2*y'+x-mu1*(x+mu)/r1^3-mu*(x-mu1)/r2^3
• y''=-2*x'+y-mu1*y/r1^3-mu*y/r2^3

?其中，mu=1/82.45,mu1=1-mu,r1=sqrt((x+mu)^2+y^2),r2=sqrt((x-mu1)^2+y^2)

代碼如下：

函數程序

function f=c_7_3_fun % 定義一個接口，通過匿名函數調用 % 通過定義接口的方式，在一個函數文件例里可編寫多個函數 f.apolloeq=@apolloeq; end % 選擇一組狀態向量：x1=x,x2=x',x3=y,x4=y' function dx=apolloeq(t,x) % 帶有附加參數的微分方程描述曲線 mu=1/82.45; mu1=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mu1)^2+x(3)^2); dx=[x(2);2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3;x(4);-2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3]; end

主程序

%% 高階常微分方程組變換方法 f=c_7_3_funx0=[1.2;0;0;-1.04935751]; tic; [t,y]=ode45(@f.apolloeq,[0,20],x0); toc; length(t) % 讀取向量長度 plot(y(:,1),y(:,3)); % 繪制相平面圖,默認精度不夠（0.001）%精度改為1e-5 tic; [t1,y1]=ode45(@f.apolloeq,[0,20],x0,odeset('RelTol',1e-5)); toc; length(t1) % 讀取向量長度 figure; plot(y1(:,1),y1(:,3)); %精度改為1e-6 tic; [t2,y2]=ode45(@f.apolloeq,[0,20],x0,odeset('RelTol',1e-6)); toc; length(t2) % 讀取向量長度 figure; plot(y2(:,1),y2(:,3)); min(diff(t2)) %計算步長 figure; plot(t2(1:end-1),diff(t2)); % 繪制實際使用的步長變化曲線 %% 采用定步長四階RK算法 f=c_7_3_funx0=[1.2;0;0;-1.04935751]; tic; [t3,y3]=rk_4(@f.apolloeq,[0,20,0.01],x0); toc; length(t3) % 讀取向量長度figure; plot(y3(:,1),y3(:,3)); % 步長改為0.001 f=c_7_3_funx0=[1.2;0;0;-1.04935751]; tic; [t4,y4]=rk_4(@f.apolloeq,[0,20,0.001],x0);toc;length(t4) % 讀取向量長度figure; plot(y4(:,1),y4(:,3)); %在求解常微分方程時建議采用變步長算法

運行結果：（算法比較）

算法：ode45()，默認精度0.001，歷時 0.007719 秒，ans =689（向量長度）

算法：ode45()，精度1e-5，歷時 0.017587 秒，ans =1309（向量長度）

算法：ode45()，精度1e-6，歷時 0.017856 秒，ans =1873（向量長度），ans =1.8927e-04（最小步長）

?注：改變算法精度重新求解，發現結果是否不同，可以通過此法驗證數值解的正確性。

算法：定步長四階RK算法，步長0.01，歷時 0.044436 秒，ans =2001（向量長度）結果錯誤

?算法：定步長四階RK算法，步長0.001，歷時 0.716064 秒，ans =20001（向量長度）

注：定步長四階RK算法中，當步長減小至0.001時，結果才大致正確，大求解時間明顯變長，建議采用變步長算法。

例1.3?將下列方程組轉化為一階微分方程組。

• x''+2*y'*x=2*y''
• x''*y'+3*x'*y''+x*y'-y=5?

代碼如下：

%% 高階微分方程隱式地含有最高項 % 除了手工轉化，可以利用MATLAB符號運算工具箱 syms x1 x2 x3 x4 dx dy % 聲明符號變量 [dx,dy]=solve([dx+2*x4*x1==2*dy,dx*x4+3*x2*dy+x1*x4-x3==5],[dx,dy]) % 解代數方程 % solve(eqn,var)

運行結果：

dx =(2*(x3 - x1*x4 - 3*x1*x2*x4 + 5))/(3*x2 + 2*x4)
dy =(2*x1*x4^2 - x1*x4 + x3 + 5)/(3*x2 + 2*x4)

注：可以選擇手動求解，正常選擇狀態變量，然后代入法求x''和y''，也可利用solve求代數方程的方法，更高效。?

剛性微分方程求解

剛性微分方程：一些解變化緩慢，另一些變化快，且相差懸殊。

注：剛性問題——在用微分方程描述的一個變化過程中，往往又包含著多個相互作用但變化速相差十分懸殊的子過程，這樣一類過程就認為具有“剛性”。描述這類過程的微分方程初值問題稱為“剛性問題”。例如，宇航飛行器自動控制系統一般包含兩個相互作用但效應速度相差十分懸殊的子系統，一個是控制飛行器質心運動的系統，當飛行器速度較大時，質心運動慣性較大，因而相對來說變化緩慢；另一個是控制飛行器運動姿態的系統，由于慣性小，相對來說變化很快，因而整個系統就是一個剛性系統。

例2.1?（例1.1）求y''+mu*(y^2-1)*y'+y=0，初值為y(0)=0.2,y'(0)=-0.7,用剛性方程解法求解mu=1000，tn=3000，x0=[2;0]的情況。

常用求解器：ode15s()

代碼如下：

%% 剛性微分方程求解1% 對于剛性問題的求解，ode15s求解速度比ode45效率高得多ff=odeset;ff.RelTol=1e-10;ff.AbsTol=1e-10; % 設置控制精度x0=[2;0];tn=3000;f=@(t,x,mu) [x(2);-mu*(x(1)^2-1)*x(2)-x(1)];mu=1000;tic;[t4,y4]=ode15s(f,[0,tn],x0,ff,mu);toc;length(t4)plot(t4,y4(:,1));figure;plot(t4,y4(:,2),'--');

運行結果：（mu=1000，tn=3000，x0=[2;0]時）

算法：ode15s，實線為[t,y]；虛線為[t,y']；用時為0.064681 秒，取點5917個。

對比組（代碼見例1.1）

算法：ode45，實線為[t,y]；虛線為[t,y']；用時為28.018461 秒，取點6757169個。

?注：1、可見，在對剛性問題的求解中，用專用的剛性求解器（此處為ode15s）能大大提高效率。2、并不是所有剛性問題求解都要刻意選擇鋼性問題的解法，許多傳統的剛性問題采用普通求解方法就可以解出。

例2.2?考慮下列常微分方程

• y'1=0.04*(1-y1)-(1-y1)*y1+0.0001*(1-y2)^2
• y'2=-10^4*y1+3000*(1-y2)^2

其中，初值為y1(0)=0,y2(0)=1，計算區間為(0,100)，選擇適當的方法求數值解。

代碼如下：

%% 剛性微分方程求解2% 用ode45求解f=@(t,y)[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;... -10^4*y(1)+3000*(1-y(2))^2]; % 用匿名函數描述微分微分方程tic;[t2,y2]=ode45(f,[0,100],[0;1]);toc;length(t2)figure;plot(t2,y2)% 改用ode15s求解f=@(t,y)[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;... -10^4*y(1)+3000*(1-y(2))^2]; % 用匿名函數描述微分微分方程tic;[t2,y2]=ode15s(f,[0,100],[0;1]);toc;length(t2)figure;plot(t2,y2)% 改用ode15s求解,提高精度ff=odeset;ff.RelTol=1e-10;ff.AbsTol=1e-10; % 設置控制精度f=@(t,y)[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;... -10^4*y(1)+3000*(1-y(2))^2]; % 用匿名函數描述微分微分方程tic;[t2,y2]=ode15s(f,[0,100],[0;1],ff);toc;length(t2)figure;plot(t2,y2)

運行結果：

分別采用ode45、ode15s、提高精度后的ode15s求解。

ode45（默認精度0.001）：歷時 0.486060 秒，ans =?356941（點數）

ode15s（默認精度0.001）：歷時 0.003206 秒，ans =56（點數）

ode15s（1e-10）：歷時 0.012131 秒，ans =614（點數）

隱式微分方程求解

隱式微分方程：不能轉換成一階顯式微分方程的微分方程。

常用求解器：ode15i()，ode15s()（代數方程）

解法：

1、矩陣法

例3.1?給定隱式微分方程組，已知x1(0)=x2(0)=0，求數值解。

• x'1*sin(x1)+x'2*cos(x2)+x1=1
• -x'1*cos(x2)+x'2*sin(x1)+x2=0

解：將其改寫成矩陣形式A(x)*x'=B(x)

• A(x)=[sin(x(1)) cos(x(2));-cos(x(2)) sin(x(1))]
• B(x)=[1-x(1);-x(2)]

假設A(x)為非奇異矩陣，則x'=A(X)^(-1)*B(x)（該方法只適合A(x)為非奇異矩陣的情況）

代碼如下：

%% 隱式微分方程求解% 若得不到奇異矩陣的報錯，則正確f=@(t,x)inv([sin(x(1)) cos(x(2));-cos(x(2)) sin(x(1))])*[1-x(1);-x(2)];% inv矩陣求逆opt=odeset;opt.RelTol=1e-6;[t,x]=ode45(f,[0,10],[0;0],opt);plot(t,x)% 此解法有隱患，參照代數方程解法可以更好的求解

運行結果：

注：此解法有隱患（需要非奇異矩陣的條件），參照下一章 “代數方程解法”?可以更好的求解（例4.2）

2、巧用代數方程

例3.2?考慮下列復雜隱式方程

• x''sin(y')+y''^2=-2*x*y*e^(-x')+x*x''*y'
• x*x''*y''+cos(y'')=3*y*x'*e^(-x)

選定狀態變量x1=x,x2=x',x3=y,x4=y'，設p1=x'',p2=y''，原方程組變為：

• p(1)*sin(x(4))+p(2)^2+2*x(1)*x(3)*exp(-x(2))-x(1)*p(1)*x(4)=0
• x(1)*p(1)*p(2)+cos(p(2))-3*x(3)*x(2)*exp(-x(1)) =0

?代碼如下：

函數程序

function f=c_7_3_fun% 定義一個接口，通過匿名函數調用% 通過定義接口的方式，在一個函數文件例里可編寫多個函數%f.apolloeq=@apolloeq;f.c7impode=@c7impode;endfunction dy=c7impode(t,x) % x為參數dx=@(p,x) [p(1)*sin(x(4))+p(2)^2+2*x(1)*x(3)*exp(-x(2))-x(1)*... p(1)*x(4);x(1)*p(1)*p(2)+cos(p(2))-3*x(3)*x(2)*exp(-x(1))]; % 用x1,x2,x3,x4表示p1,p2ff=optimset;ff.Display='off';% Level of display (see Iterative Display): % 'off' or 'none' displays no output. % 'iter' displays output at each iteration, % and gives the default exit message. % 'iter-detailed' displays output at each iteration, % and gives the technical exit message. % 'final' (default) displays just the final output, % and gives the default exit message. % 'final-detailed' displays just the final output, % and gives the technical exit message.dx1=fsolve(dx,x([1,3]),ff,x); % x = fsolve(fun,x0,options) % fsolve求解非線性方程，初值選擇p1(0)=x1,p2(0)=x3,加快代數方程收斂速度和精度 % x作為已知附加參數給出，得出結果為p1,p2dy=[x(2);dx1(1);x(4);dx1(2)]; % dy=[x2=x',dx1(1)=p1=x'',x4=y',dx1(2)=p2=y''] % 求ode對應的解為[x,x',y,y']end

主程序

%% 復雜隱式微分方程f=c_7_3_fun[t,x]=ode15s(f.c7impode,[0,2],[1,0,0,1]);% x1隨t變化的圖像figure(1);plot(t,x(:,1));% x2隨t變化的圖像figure(2);plot(t,x(:,2));% x3隨t變化的圖像figure(3);plot(t,x(:,3));% x4隨t變化的圖像figure(4);plot(t,x(:,4));% x隨t變化的總圖像figure(5);plot(t,x);% 下面代碼找不到顯式解（方程和系統求解器無解）% syms x1 x2 x3 x4 dx dy % 聲明符號變量% [dx,dy]=solve([dx*sin(x4)+dy^2+2*x1*x3*exp(-x2)-x1*...% dx*x4==0,x1*dx*dy+cos(dy)-3*x3*x2*exp(-x1)==0],[dx,dy])

?運行結果：

?

3、直接用ode15i()函數求解

?求解之前需要給出（x0,x'0），不能隨意賦值，只能有n個是獨立的，其余的需要用隱式方程求解，否則將出現矛盾的初始條件，用decic()函數能的處相容的初值。

例 3.3 （例3.2）考慮下列復雜隱式方程

• x''sin(y')+y''^2=-2*x*y*e^(-x')+x*x''*y'
• x*x''*y''+cos(y'')=3*y*x'*e^(-x)

選定狀態變量x1=x,x2=x',x3=y,x4=y'，原方程組變為隱式方程標準型：

• x'1-x2=0
• x'2*sin(x4)+x'4^2+2*e^(-x2)*x1*x3-x1*x'2*x4=0
• x'3-x4=0
• x1*x'2*x'4+cosx'4-3*e^(-x1)*x3*x2=0

代碼如下：

%% 利用ode15i()求解復雜隱式微分方程% 隱式微分方程的描述f=@(t,x,xd)[xd(1)-x(2); xd(2)*sin(x(4))+xd(4)^2+2*exp(-x(2))*x(1)*x(3)-x(1)*xd(2)*x(4); xd(3)-x(4); x(1)*xd(2)*xd(4)+cos(xd(4))-3*exp(-x(1))*x(3)*x(2)];x0=[1,0,0,1]; % 初值xd0=[1,2,3,-5]; % 無需給出確定值，但要判斷正負x0f=[1 1 1 1]; % 保留x0xd0f=[]; % 重新計算xd[x0,xd0]=decic(f,0,x0,x0f,xd0,xd0f)% decic為ode15i計算一致的初始條件% t0=0% 求出相容的初值,由x0確定x0'[t,x]=ode15i(f,[0,2],x0,xd0);plot(t,x) % 繪制時間響應曲線

?運行結果：（同例3.2）

x0 =? ? ? ?1? ? ?0? ? ?0? ? ?1
xd0 =? -0.0000? 1.6833? ? 1.0000? ?-0.5166

微分代數方程求解

微分代數方程：指在微分方程中，某些變量間滿足滿足某些代數方程的約束。

常用求解器：ode15s()、ode45()

微分方程更一般的形式可以寫成

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

對于真正的微分代數方程，M(t,x)為奇異矩陣

例4.1?考慮下面微分代數方程

• x'1=-0.2*x1+x2*x3+0.3*x1*x2
• x'2=2*x1*x2-5*x2*x3-2*x2*x2
• 0=x1+x2+x3-1

?初始條件：x1(0)=0.8,x2(0)=x3(0)=0.1

解：用矩陣形式表示該微分代數方程

[1 0 0;0 1 0;0 0 0]*[x'1 x'2 x'3]=[-0.2*x1+x2*x3+0.3*x1*x2;2*x1*x2-5*x2*x3-2*x2*x2;x1+x2+x3-1]

代碼如下：

%% 微分代數方程求解f=@(t,x)[-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2); 2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2); x(1)+x(2)+x(3)-1];M=[1,0,0;0,1,0;0,0,0]; % 輸入質量矩陣options=odeset;options.Mass=M; % 設置質量矩陣x0=[0.8;0.1;0.1];[t,x]=ode15s(f,[0,20],x0,options);plot(t,x)% 也可以轉換成常微分方程求解fDae=@(t,x)[-0.2*x(1)+x(2)*(1-x(1)-x(2))+0.3*x(1)*x(2); 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)]; % 方程右側的形式模型x0=[0.8,0.1]; % 初值[t,x]=ode45(fDae,[0,20],x0);x3=1-x(:,1)-x(:,2);figure;plot(t,x,t,x3);% 也可以用隱式微分方程求解的方法求解f=@(t,x,xd)[xd(1)+0.2*x(1)-x(2)*x(3)-0.3*x(1)*x(2); xd(2)-2*x(1)*x(2)+5*x(2)*x(3)+2*x(2)*x(2); x(1)+x(2)+x(3)-1]; % 隱式微分方程x0=[0.8,0.1,2];x0f=[1,1,0];xd0=[1,1,1];xd0f=[];[x0,xd0]=decic(f,0,x0,x0f,xd0,xd0f); % 相容初始條件% 隱式微分方程求解與繪圖res=ode15i(f,[0,20],x0,xd0);% res為struct結構體數組plot(res.x,res.y)

運行結果：

例4.2 （例3.1）?給定隱式微分方程組，已知x1(0)=x2(0)=0，求數值解。

• x'1*sin(x1)+x'2*cos(x2)+x1=1
• -x'1*cos(x2)+x'2*sin(x1)+x2=0

?代碼如下：

f=@(t,x)[1-x(1);-x(2)];M=@(t,x)[sin(x(1)),cos(x(2));-cos(x(2)),sin(x(1))]; % 質量矩陣options=odeset;options.Mass=M;options.RelTol=1e-6; % 設置求解精度[t,x]=ode45(f,[0,10],[0;0],options);plot(t,x)

運行結果：

總結

以上是生活随笔為你收集整理的几类常微分方程的matlab求解方法 | 刚性微分方程、隐式微分方程、微分代数方程的全部內容，希望文章能夠幫你解決所遇到的問題。

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