MATLAB解常微分方程組的解法(好東西要共享)


1:問題

常微分方程的初值問題的標准數學表述為:y'=f(t,y),a<=t<=b,y(a)=y(0)

 

;我們要求解的任何高階常微分方程都可以用替換法化為上式所示的一階形式,其中y為向量,yo為初始值。

2:Matlab中解決以上問題的步驟

   (1):化方程組為標准形式。

          例如:y'''-3y''-y’y=0,y(0)=0,y'(0)=1,y''(0)=-1.

          把微分方程的高階導數寫為低階導數的算式,即:

           y'''=3y''+y'y,設:y1=y,y2=y',y3=y'',則原方程化為下列等價的方程組:

          滿足初值條件: 已把該方程化成了標准形式。

其中:y'->(y1’,y2’,y3’),a->(0,0,0),y0->(0,1,-1),f(t,y)->(y2,y3,3y3+y2y1).

(2):把微分方程組編成m函數文件。

      如:function dy=F(t,y)

dy=[y(2);y(3);3*y(3)+y(2)*y(1)];

注意:A:在函數文件里,雖然寫微分方程時並不同時包含參數t和y,但第一行必須包含這兩個輸入變量。B:向量dy必須為列向量。

(3):調用一個微分方程的求解函數求解。

    [T,Y]=solver(‘F’,tspan,y0);

其中:solver:求解函數名;

F:包含微分方程的m文件; 

tspan為積分的數據范圍,其格式為:[t0,tfinal]; 

y0為t0時刻的初值列向量。

輸出參數T和Y為列向量

T為時刻向量。

Y表是不同時刻的函數值。

3:一個求解常微分方程初值問題的完整過程。

問題:求解方程y’’-3(1-y^2)y’+y=0在初值y’(0)=3,y(0)=2的解。

1化成標准形式:

設y1=y,y2=y’,則: 初值為:

2編寫函數文件ode.m,類容為:

function dy=ode(t,y)

dy=[y(2);3*(1-y(1)^2)*y(2)-y(1)];

3調用函數ode45求解,時間區間為[0,20]:

[T,Y]=ode45(‘ode’,[0,20],[2;3]);

輸出結果[T,Y]中T為時間點組成的向量。Y為對應於T中時間點的y(1)和y(2)的值。

4繪制解的曲線,結果如圖。

plot(T,Y(:,1),’-’,T,Y(:,2),’--’)

title(‘Solution of ODE Equation’);

xlabel(‘time T’)

ylabel(‘solution Y’);

legend(‘Y1’,’Y2’)

Matlab利用數值方法來求解常微分方程的解,其思路如下:把求解的時間區間划分成有限步,對應於每一步將計算出一個解,如果求得的解不滿足誤差限制,則減少步長,再求解。如此重復,直到滿足誤差限為止。

a剛性問題(stiff):方程組的解不同分量的數量級差別較大,對於數值求解是一大困難。Matlab既能解決非剛性問題,也能解決剛性問題。

b三個解決非剛性問題的函數:ode45,ode23,ode113

c兩個解剛性問題的函數:ode15s和ode23s


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM