當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

四阶龙格库塔法(Runge-Kutta)求解常微分方程的 Matlab程序及案例

發布時間：2023/12/31 编程问答 37 豆豆

生活随笔收集整理的這篇文章主要介紹了四阶龙格库塔法(Runge-Kutta)求解常微分方程的 Matlab程序及案例小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

文章目錄

- 1. 算法
- 2. 程序
- 3. 案例
- 4. 聯系作者

1. 算法

上一篇介紹了顯式歐拉法、隱式歐拉法、兩步歐拉法和改進歐拉法求解常微分方程初值問題；其中顯式歐拉法和隱式歐拉法是一階算法精度，截斷誤差為 $O(h2)O\left( {{h^2}} \right)$ ；兩步歐拉法和改進歐拉法是二階算法精度，截斷誤差為 $O(h3)O\left( {{h^3}} \right)$ ；歐拉法的精度有限、需要求解步長 $h$ 很小。本篇介紹求解精度更高的四階龍格庫塔法(Runge-Kutta)，其截斷誤差為 $O(h5)O\left( {{h^5}} \right)$ 。

對于一階微分方程初值問題：

${y˙=f(t,y)y(t0)=y0\left\{ \begin{array}{l} {\bf{\dot y}} = f\left( {t,{\bf{y}}} \right) \\ {\bf{y}}\left( {{t_0}} \right) = {{\bf{y}}_0} \\ \end{array} \right.$

式中， ${t_0}$ 為初始時間（已知常數）， $y0{{\bf{y}}_0}$ 為初始狀態（已知向量）， $f(t,y)f\left( {t,{\bf{y}}} \right)$ 為關于時間 ${t}$ 和狀態 $y{{\bf{y}}}$ 的函數（已知函數）。

四階龍格庫塔法(Runge-Kutta)求解算法為：

$k1=f(t(k),y(k)){{k}_{1}}=f\left( t\left( k \right),\mathbf{y}\left( k \right) \right)$
$k2=f(t(k)+h2,y(k)+h2k1){{k}_{2}}=f\left( t\left( k \right)+\frac{h}{2},\mathbf{y}\left( k \right)+\frac{h}{2}{{k}_{1}} \right)$
$k3=f(t(k)+h2,y(k)+h2k2){{k}_{3}}=f\left( t\left( k \right)+\frac{h}{2},\mathbf{y}\left( k \right)+\frac{h}{2}{{k}_{2}} \right)$
$k4=f(t(k)+h,y(k)+hk3){{k}_{4}}=f\left( t\left( k \right)+h,\mathbf{y}\left( k \right)+h{{k}_{3}} \right)$
$y(k+1)=y(k)+h6(k1+2k2+2k3+k4)\mathbf{y}\left( k+1 \right)=\mathbf{y}\left( k \right)+\frac{h}{6}\left( {{k}_{1}}+2{{k}_{2}}+2{{k}_{3}}+{{k}_{4}} \right)$
$y(0)=y0\mathbf{y}\left( 0 \right)={{\mathbf{y}}_{0}}$

上式是關于 $y(k){\bf{y}}\left( k \right)$ 向 $y(k+1){\bf{y}}\left( k+1 \right)$ 的遞推形式，可以根據初始條件按照遞推關系依次求出 $y(1),y(2),y(3),y(4)?,y(N)?{\bf{y}}\left( 1 \right),{\bf{y}}\left( 2 \right),{\bf{y}}\left( 3 \right),{\bf{y}}\left( 4 \right) \cdots ,{\bf{y}}\left( N \right) \cdots$ ，此離散序列即為微分方程的數值解。

微分方程的數值解法，本質是使用數值積分來實現 $y˙{\bf{\dot y}}$ 向 $y{\bf{y}}$ 的轉換。四階龍格庫塔法通過對微分的四步分段逼近，在一個求解步長內能夠逼近復雜的曲線，因此能夠取得較高的計算精度，其截斷誤差為 $O(h5)O\left( {{h^5}} \right)$ 。

2. 程序

作者使用Matlab開發了四階龍格庫塔法求解常微分方程的程序，能夠方便快捷的求解一階常微分方程的初值問題。

function [T,X,dX] = ODE_RK4( Hfun,t,h,x0 ) % [T,X] = ODE_RK4( Hfun,t,h,x0 ) 4階龍格-庫塔法求解常微分方程 % Hfun為描述狀態導數的函數句柄，格式為 dX = Hfun( t,X ) % 必須保證返回dX的格式為行向量 % t為時間節點，可為標量，時間范圍為 T = 0:h:t % 長2向量，時間范圍為 T = t(1):h:t(2) % 向量，時間范圍為 T = t % h為時間步長，在t的前兩種情況下，必須給出h具體值 % 直接給出時間節點t時，h可用[]或任意值占位 % x0為狀態量初始值 % 算法： % K1 = Hfun( t(k-1),X(k-1) ) = dX(k-1) % K2 = Hfun( t(k-1)+h/2,X(k-1)+h*K1/2 ) % K3 = Hfun( t(k-1)+h/2,X(k-1)+h*K2/2 ) % K4 = Hfun( t(k-1)+h ,X(k-1)+h*K3 ) % X(k) = X(k-1) + (h/6) * (K1 + 2*K2 + 2*K3 +K4) % By ZFS@wust 2021 % 獲取更多Matlab/Simulink原創資料和程序，清關注微信公眾號：Matlab Fans

下面結合實例進行演示和分析。

3. 案例

求解一階常微分方程（式中向量 $x{\bf{x}}$ 等價于前文中的向量 $y{\bf{y}}$ ）：
$x˙=f(t,x)=[x(2)?x(3)?x(1)?x(3)?0.51?x(1)?x(2)]\mathbf{\dot{x}}=f\left( t,\mathbf{x} \right)=\left[ \begin{matrix} \mathbf{x}(2)*\mathbf{x}(3) \\ -\mathbf{x}(1)*\mathbf{x}(3) \\ -0.51*\mathbf{x}(1)*\mathbf{x}(2) \\ \end{matrix} \right]$
$x(0)=[011]\mathbf{x}\left( 0 \right)=\left[ \begin{matrix} 0 \\ 1 \\ 1 \\ \end{matrix} \right]$

% 四階龍格庫塔算法(Runge-Kutta)測試程序 % By ZFS@wust 2021 % 獲取更多Matlab/Simulink原創資料和程序，清關注微信公眾號：Matlab Fansclear clc close all%% 仿真步長 h=0.6 時 Hfun = @(t,x) [ x(2)*x(3); -x(1)*x(3); -0.51*x(1)*x(2)]; % 一階微分方程導數表達式% 參數 t = [0 12]; % 時間范圍 h = 0.6; % 時間步長 x0 = [0 1 1]; % 初始狀態% 改進歐拉法求解 [T1,X1] = ODE_ImprovedEuler( Hfun,t,h,x0 ); % 改進歐拉法求解 [T2,X2] = ODE_RK4( Hfun,t,h,x0 ); % Matlab自帶ode45求解 [T3,X3] = ode45( Hfun,t(1):h:t(2),x0 );% 繪圖對比 figure subplot(311) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_1') legend('改進歐拉法','四階龍格庫塔法','ode45') subplot(312) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_2') legend('改進歐拉法','四階龍格庫塔法','ode45') subplot(313) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_3') legend('改進歐拉法','四階龍格庫塔法','ode45')%% 仿真步長 h=0.9 時 % 參數 h = 0.9; % 時間步長% 改進歐拉法求解 [T1,X1] = ODE_ImprovedEuler( Hfun,t,h,x0 ); % 改進歐拉法求解 [T2,X2] = ODE_RK4( Hfun,t,h,x0 ); % Matlab自帶ode45求解 [T3,X3] = ode45( Hfun,t(1):h:t(2),x0 );% 繪圖對比 figure subplot(311) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_1') legend('改進歐拉法','四階龍格庫塔法','ode45') subplot(312) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_2') legend('改進歐拉法','四階龍格庫塔法','ode45') subplot(313) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_3') legend('改進歐拉法','四階龍格庫塔法','ode45')

不同時間步長 $h$ 時的數值計算結果：

步長 $h = 0.6 s$
步長 $h = 0.9 s$

從計算結果可以看出，四階龍格庫塔法(Runge-Kutta)即使在步長很大時，也能保持較高的求解精度，求解精度與Matlab自帶的ode45函數相當，相對于改進歐拉算法求解精度有明顯提高。
自己編程實現四階龍格庫塔法(Runge-Kutta)，相對于直接調用ode45等Matlab自帶的龍格庫塔法的最大優勢在于：可以將求解程序和模型描述文件融合起來，解決各類參數時變、條件判斷、多模型切換等問題。后續補充實際案例。

4. 聯系作者

有Matlab/Simulink方面的技術問題，歡迎發送郵件至944077462@qq.com討論。
更多Matlab/Simulink原創資料，歡迎關注微信公眾號：Matlab Fans

源程序下載:
1. csdn資源: 四階龍格庫塔法(Runge-Kutta)求解常微分方程的Matlab程序及案例
2. 掃碼關注微信公眾號Matlab Fans，回復BK09獲取百度網盤下載鏈接。

總結

以上是生活随笔為你收集整理的四阶龙格库塔法(Runge-Kutta)求解常微分方程的 Matlab程序及案例的全部內容，希望文章能夠幫你解決所遇到的問題。

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

上一篇：帆软报表嵌入python程序_C#教程之
下一篇：双线性的定义以及他的性质