【正文】 )/(y_rk4(k)+h*k2/2)。 k4=(y_rk4(k)+h*k3)2*(t(k)+h)/(y_rk4(k)+h*k3)。 y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6。enddisp(39。 時間 解析解 歐拉法 RK2法 RK4法39。)yt=[t39。, y39。, y_euler39。, y_rk239。, y_rk439。]。disp(yt)【運行結果】 時間 解析解 歐拉法 RK2法 RK4法 0 %H= 時間 解析解 歐拉法 RK2法 RK4法 0 H= 時間 解析解 歐拉法 RK2法 RK4法 0 %歐拉法需要步長足夠小時才逼近解析解，RK4逼近得最快，其次是RK22. 考慮如下二階系統： 在上的數字仿真解（已知：，），并將不同步長下的仿真結果與解析解進行精度比較。 說明：已知該微分方程的解析解分別為： 采用RK4法進行計算，選擇狀態(tài)變量： 則有如下狀態(tài)空間模型及初值條件 采用RK4法進行計算。程序如下：clearh=input(39。請輸入步長h=39。)。 % 輸入步長M=round(10/h)。 % 置總計算步數t(1)=0。 % 置自變量初值y_0(1)=100。 y_05(1)=100。 % 置解析解的初始值(y_0和y_05分別對應于為R=0和R=)x1(1)=100。 x2(1)=0。 % 置狀態(tài)向量初值y_rk4_0(1)=x1(1)。 y_rk4_05(1)=x1(1)。 % 置數值解的初值 % 求解析解for k=1:M t(k+1)=t(k)+h。 y_0(k+1)=100*cos(t(k+1))。 y_05(k+1)=100*exp(t(k+1)/2).*cos(sqrt(3)/2*t(k+1))+100*sqrt(3)/3*exp(t(k+1)/2).*sin(sqrt(3)/2*t(k+1))。end% 利用RK4法求解% R=0for k=1:M k11=x2(k)。 k12=x1(k)。 k21=x2(k)+h*k12/2。 k22=(x1(k)+h*k11/2)。 k31=x2(k)+h*k22/2。 k32=(x1(k)+h*k21/2)。 k41=x2(k)+h*k32。 k42=(x1(k)+h*k31)。 x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6。 x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6。 y_rk4_0(k+1)=x1(k+1)。end% R=for k=1:M k11=x2(k)。 k12=x1(k)x2(k)。 k21=x2(k)+h*k12/2。 k22=(x1(k)+h*k11/2)(x2(k)+h*k12/2)。 k31=x2(k)+h*k22/2。 k32=(x1(k)+h*k21/2)(x2(k)+h*k22/2)。 k41=x2(k)+h*k32。 k42=(x1(k)+h*k31)(x2(k)+h*k32)。 x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6。 x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6。 y_rk4_05(k+1)=x1(k+1)。end% 求出誤差最大值err_0=max(abs(y_0y_rk4_0))。err_05=max(abs(y_05y_rk4_05))。% 輸出結果disp(39。最大誤差（R=0) 最大誤差（R=)39。)err_max=[err_0,err_05]。disp(err_max) ，并且將數值解直接與解析解比較，求出最大誤差數據列入下表中。步長hR=0最大誤差101010101010109108105104101R=最大誤差1011101210121010109106105102從上表中可以看出，當步長h=，總誤差最??；，由于舍入誤差變大而使總誤差增加；，則由于截斷誤差的增加也使得總誤差加大。另外，當系統的解變化激烈時（如R=0），誤差對步長的變化較為敏感；當系統的解變化平穩(wěn)時，步長的變化對誤差的影響就要緩和得多。數值積分算法確定以后，在選擇步長時，需要綜合考慮。 41