此處參考教材為李榮華的《微分方程數值解法》
使用工具：Matlab

1. 算法

注: 最后一行應為k₄,上面為筆誤

2. 算法

I.需要求解的函數

function f=f1D(t,u,ft)% ft為方程編號，u1D為精確解函數u(t),注意與f1D對應右端項函數f(t,u(t))switch ftcase 1 %P10 習題1%v=exp(-5*t);f=-5*u;case 2 % P27 例3.1%v=(1+t.^2).^2;f=4*t.*u.^(1/2);case 3 %補充練習 P54 例1%v=exp(-t);f=-u; end

II.對應精確解

function v=u1D(t,ft)% ft為方程編號，u1D為精確解函數u（t）,注意與f1D 對應右端項函數f(t,u(t))switch ftcase 1 %P10 習題1v=exp(-5*t);%f=-5*u;case 2 % P27 例3.1v=(1+t.^2).^2;%f=4*t.*u.^(1/2);case 3 %補充練習 P54 例1%v=exp(-t);f=-u; end

III.求解問題的主函數

function ODESolver(ft,t0,T,h)t=(t0:h:T);%t節點 n=length(t);%步數 ExactU=u1D(t,ft);%精確解 %-------------------------------------- % 1)向前歐拉，單步顯式，u(i)=u(i-1)+h*f(t(i-1),u(i-1)) u=0.*t;%數值解初始化 u(1)=ExactU(1);%設初值 tic; %記錄CPU計算開始時間 for i=2:nu(i)=u(i-1)+h*f1D(t(i-1),u(i-1),ft); end time1=toc; %記錄CPU計算結束時間 NumU1=u;%計算結束賦值 Err1=(NumU1(n)-ExactU(n));%T點處誤差 RelErr1=Err1/ExactU(n);%相對誤差%-------------------------------------- % 2)向后Euler法，單步隱式，u(i)=u(i-1)+h*f(t(i),u(i)) u=0.*t;%數值解初始化 u(1)=ExactU(1);%設初值 tic; %記錄CPU計算開始時間 for i=2:nuc1=u(i-1)-1;%非線性迭代設初值uc2=u(i-1);%非線性迭代設初值，因為迭代uc1=uc2,故初值uc2=u(i-1)delta=10^(-5);%迭代閾值while abs(uc2-uc1)>deltauc1=uc2;uc2=u(i-1)+h*f1D(t(i),uc1,ft);endu(i)=uc2;%非線性迭代結束 end time2=toc; %記錄CPU計算結束時間 NumU2=u; %計算結束賦值 Err2=NumU2(n)-ExactU(n); %T點處誤差 RelErr2=Err2/ExactU(n); %相對誤差%-------------------------------------- % 3）梯形法，單步隱式法， u(i)=u(i-1)+(h/2)*(f(t(i-1),u(i-1))+f(t(i),u(i))) u=0.*t; %數值解初始化 u(1)=ExactU(1); %設初值 tic; %記錄CPU計算開始時間 for i=2:nuc1=u(i-1)-1; %非線性迭代設置初值uc2=u(i-1); %非線性迭代設置初值，因為迭代uc1=uc2,故初值uc2=u(i-1)delta=10^(-5); %d迭代閾值while abs(uc2-uc1)>deltauc1=uc2;uc2=u(i-1)+(h/2)*(f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));endu(i)=uc2; %非線性迭代結束 end time3=toc; %記錄CPU計算結束時間 NumU3=u; %計算結束賦值 Err3=NumU3(n)-ExactU(n); %T點處誤差 RelErr3=Err3/ExactU(n); %相對誤差%-------------------------------------- % 4) 改進Euler法，單步顯示法，uc（i）=u(i-1)+h*f(t(i-1),u(i-1)),u(i)=u(i-1)+(h/2)*(f(t(i-1),u(i-1))+f(t(i),uc(i))) u=0.*t;%數值解法初始化 u(1)=ExactU(1);%設初值 tic;%記錄CPU計算開始時間 for i=2:nuc1=u(i-1)+h*4*t(i-1)*(u(i-1)^(1/2));%用向前Euler法預估u(i)=u(i-1)+(h/2)*(f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));%用梯形法校正 end time4=toc;%記錄CPU計算結束時間 NumU4=real(u);%計算結束賦值 Err4=NumU4(n)-ExactU(n);%T點處誤差 RelErr4=Err4/ExactU(n);%相對誤差%-------------------------------------- %5）二步法a=0,二步顯式法,u(n+2)-(1+a)*u(n+1)+a*u(n)=(1/2)*h*((3-a)*f(n+1)-(1+a)*f(n)) u=0.*t;%數值解初始化 u(1)=ExactU(1);%設初值 u(2)=ExactU(2);%設初值 a=0; tic;%記錄CPU計算開始時間 for i=3:n u(i)=(1+a)*u(i-1)-a*u(i-2)+(1/2)*h*((3-a)*f1D(t(i-1),u(i-1),ft)-(1+a)*f1D(t(i-2),u(i-2),ft)); end time5=toc;%記錄CPU計算結束時間 NumU5=u;%計算結束賦值 Err5=NumU5(n)-ExactU(n);%T點處誤差 RelErr5=Err5/ExactU(n);%相對誤差%-------------------------------------- %6)二步法a=-5,二步顯式法，u(n+2)-(1+a)*u(n+1)+a*u(n)=(1/2)*h*((3-a)*f(n+1)-(1+a)*f(n)) u=0.*t;%數值解初始化 u(1)=ExactU(1); %設初值 u(2)=ExactU(2); %設初值 a=-5; tic; %記錄CPU計算結束時間 for i=3:n u(i)=(1+a)*u(i-1)-a*u(i-2)+(1/2)*h* ((3-a)*f1D(t(i-1),u(i-1),ft)-(1+a)*f1D(t(i-2),u(i-2),ft)); end time6=toc; %記錄CPU計算結束時間 NumU6=real(u); %計算結束賦值，由于u計算出現負數，開根號出現復數，故取實部 Err6=NumU6(n)-ExactU(n); %T點處誤差 RelErr6=Err6/ExactU(n); %相對誤差%-------------------------------------- % 7）Simpson法，二步隱式法，u(i)=u(i-2)+(2*h/6)*(f(t(i-2),u(i-2))+4*f(t(i-1),u(i-1))+f(t(i),u(i))) u=0.*t; %數值解初始化 u(1)=ExactU(1); %設初值 u(2)=ExactU(2); %設初值 tic; %記錄CPU計算結束時間 for i=3:nuc1=u(i-2); %非線性迭代設初值uc2=u(i-1); %非線性迭代設初值，因為迭代uc1=uc2，故初值uc2=u(i-1)delta=10^(-5); %迭代閾值，精度及運算時間受其影響while abs(uc2-uc1)>deltauc1=uc2; uc2=u(i-2)+(h/3)*(f1D(t(i-2),u(i-2),ft)+4*f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));endu(i)=uc2; %非線性迭代結束 end time7=toc; %記錄CPU計算結束時間 NumU7=u; %計算結束賦值 Err7=NumU7(n)-ExactU(n);%T點處誤差 RelErr7=Err7/ExactU(n);%相對誤差%-------------------------------------- % 8)四級四階RK法，單步顯式法 u=0.*t;%數值解初始化 u(1)=ExactU(1);%設初值 tic;%記錄CPU開始計算的時間 for i=2:nk1=f1D(t(i-1),u(i-1),ft);k2=f1D(t(i-1)+h/2,u(i-1)+(h/2)*k1,ft);k3=f1D(t(i-1)+h/2,u(i-1)+(h/2)*k2,ft);k4=f1D(t(i-1)+h,u(i-1)+h*k3,ft);u(i)=u(i-1)+(h/6)*(k1+2*k2+2*k3+k4); end time8=toc;%記錄CPU計算結束時間 NumU8=u;%計算結束賦值 Err8=(NumU8(n)-ExactU(n));%T點處誤差 RelErr8=Err8/ExactU(n);%相對誤差%-------------------------------------- %畫圖 figure; plot(t,ExactU,'k-*',t,NumU1,'r-o',t,NumU2,'b-o',t,NumU3,'g-o',t,NumU4,'c-o',t,NumU5,'r-d',t,NumU7,'g-d',t,NumU8,'b-d'); legend('精確解','1)向前歐拉','2)向后歐拉','3)梯形法','4)改進Euler','5)二步法a=0','7)Simpson法','8）RK法','Location','NorthEastOutside');%制表 f=figure; data=([t;ExactU;NumU1;NumU2;NumU3;NumU4;NumU5;NumU6;NumU7;NumU8])'; data=[data;[0,0,Err1,Err2,Err3,Err4,Err5,Err6,Err7,Err8];[0,0,RelErr1,RelErr2,RelErr3,RelErr4,RelErr5,RelErr6,RelErr7,RelErr8];[0,0,time1,time2,time3,time4,time5,time6,time7,time8]]; colnames={'t','精確解','1)向前歐拉','2)向后歐拉','3)梯形法','4)改進Euler','5)二步法a=0','6)二步法a=-5','7)Simpson法','8）RK法'}; tab=uitable(f,'Data',data,'ColumnName',colnames,'Position',[50 50 800 500]);

3.給出簡單結果

[^此處圖略，僅給出簡單數字結果；感興趣的讀者可以復制代碼運行一下]

P10頁：習題1

（1） 步長0.1：

在命令行窗口輸入
>>ODESolver(1,0,1,0.1)

所得最后一步誤差結果：

Err1 =-0.0058
Err2 =0.0106
Err3 =-6.9002e-04
Err4 =-0.0098
Err5 =0.0045
Err6 =-2.1480e+04
Err7 =-3.0847e-04
Err8 =2.6728e-05

（2）步長0.05：

在命令行窗口輸入
>>ODESolver(1,0,1,0.05)

所得結果：

Err1 =-0.0036
Err2 =0.0048
Err3 =-1.7339e-04
Err4 =-0.0069
Err5 = 9.8516e-04
Err6 = -2.1437e+10
Err7 =-1.2741e-05
Err8 = 1.3516e-06

P27頁：例題3.1

（1）步長0.1：

在命令行窗口輸入
>>ODESolver(2,0,5,0.1)
所得結果：

Err1 = -3.5381
Err2 = 4.0508
Err3 = 0.0806
Err4 = -0.0763
Err5 = -0.3675
Err6 = -6.9711e+08
Err7 = -1.3874e-06
Err8 = -9.5303e-05

（2）步長0.05：

在命令行窗口輸入
>>ODESolver(2,0,2,0.05)
所得結果：

Err1 = -1.8298
Err2 = 1.9580
Err3 = 0.0201
Err4 = -0.0196
Err5 = -0.0963
Err6 = -5.4646e+21
Err7 = -1.4567e-06
Err8 = -6.2349e-06

P54 例1

（1）步長0.1：

在命令行窗口輸入
>>ODESolver(3,0,2,0.1)
所得結果：

Err1 = -0.0138
Err2 = 0.0133
Err3 = -2.2421e-04
Err4 = -0.0974
Err5 = 0.0011
Err6 = -1.2434e+08
Err7 = 1.8151e-08
Err8 = 2.4519e-07

4.對算例結果進行簡單的分析

從整體來看：

第一道題在靠近1的附近，數值解精確度都比較高。第二道題在靠近0的附近，數值解精確度都比較高。其中，向前Euler法和向后Euler法落在精確解的兩側。圖表的最后三行分別表示最后一步的數值誤差、最后一步的相對誤差、每個算法所花費的時間。

從步長來看：

圖像中可以較直觀的發現，h=0.05比h=0.1的精度整體更高一些。當步長變小時，對結果的分析更加清晰。

從八種不同的算法來看：

利用直觀的圖像以及數據表，Simpson方法誤差最小，其次是RK法，二步法的誤差最大（因為這種算法不穩定，可從圖表中看出當a=-5時，二步法顯示出不穩定性，且當步長縮小，時間增加時，不穩定性表現得非常明顯）。在所有的單步法中，梯形法和改進歐拉法效果也比較好。

總結

以上是生活随笔為你收集整理的微分方程数值解法（1）——常微分方程初值问题的数值解法的全部內容，希望文章能夠幫你解決所遇到的問題。

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