當(dāng)前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

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

發(fā)布時(shí)間：2023/12/31 编程问答 44 豆豆

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

文章目錄

- 1. 算法
- 2. 程序
- 3. 案例
- 4. 聯(lián)系作者

1. 算法

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

對(duì)于一階微分方程初值問題：

${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}$ 為初始時(shí)間（已知常數(shù)）， $y0{{\bf{y}}_0}$ 為初始狀態(tài)（已知向量）， $f(t,y)f\left( {t,{\bf{y}}} \right)$ 為關(guān)于時(shí)間 ${t}$ 和狀態(tài) $y{{\bf{y}}}$ 的函數(shù)（已知函數(shù)）。

四階龍格庫塔法(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}}$

上式是關(guān)于 $y(k){\bf{y}}\left( k \right)$ 向 $y(k+1){\bf{y}}\left( k+1 \right)$ 的遞推形式，可以根據(jù)初始條件按照遞推關(guān)系依次求出 $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$ ，此離散序列即為微分方程的數(shù)值解。

微分方程的數(shù)值解法，本質(zhì)是使用數(shù)值積分來實(shí)現(xiàn) $y˙{\bf{\dot y}}$ 向 $y{\bf{y}}$ 的轉(zhuǎn)換。四階龍格庫塔法通過對(duì)微分的四步分段逼近，在一個(gè)求解步長(zhǎng)內(nèi)能夠逼近復(fù)雜的曲線，因此能夠取得較高的計(jì)算精度，其截?cái)嗾`差為 $O(h5)O\left( {{h^5}} \right)$ 。

2. 程序

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

function [T,X,dX] = ODE_RK4( Hfun,t,h,x0 ) % [T,X] = ODE_RK4( Hfun,t,h,x0 ) 4階龍格-庫塔法求解常微分方程 % Hfun為描述狀態(tài)導(dǎo)數(shù)的函數(shù)句柄，格式為 dX = Hfun( t,X ) % 必須保證返回dX的格式為行向量 % t為時(shí)間節(jié)點(diǎn)，可為標(biāo)量，時(shí)間范圍為 T = 0:h:t % 長(zhǎng)2向量，時(shí)間范圍為 T = t(1):h:t(2) % 向量，時(shí)間范圍為 T = t % h為時(shí)間步長(zhǎng)，在t的前兩種情況下，必須給出h具體值 % 直接給出時(shí)間節(jié)點(diǎn)t時(shí)，h可用[]或任意值占位 % x0為狀態(tài)量初始值 % 算法： % 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原創(chuàng)資料和程序，清關(guān)注微信公眾號(hào)：Matlab Fans

下面結(jié)合實(shí)例進(jìn)行演示和分析。

3. 案例

求解一階常微分方程（式中向量 $x{\bf{x}}$ 等價(jià)于前文中的向量 $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)測(cè)試程序 % By ZFS@wust 2021 % 獲取更多Matlab/Simulink原創(chuàng)資料和程序，清關(guān)注微信公眾號(hào)：Matlab Fansclear clc close all%% 仿真步長(zhǎng) h=0.6 時(shí) Hfun = @(t,x) [ x(2)*x(3); -x(1)*x(3); -0.51*x(1)*x(2)]; % 一階微分方程導(dǎo)數(shù)表達(dá)式% 參數(shù) t = [0 12]; % 時(shí)間范圍 h = 0.6; % 時(shí)間步長(zhǎng) x0 = [0 1 1]; % 初始狀態(tài)% 改進(jìn)歐拉法求解 [T1,X1] = ODE_ImprovedEuler( Hfun,t,h,x0 ); % 改進(jìn)歐拉法求解 [T2,X2] = ODE_RK4( Hfun,t,h,x0 ); % Matlab自帶ode45求解 [T3,X3] = ode45( Hfun,t(1):h:t(2),x0 );% 繪圖對(duì)比 figure subplot(311) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_1') legend('改進(jìn)歐拉法','四階龍格庫塔法','ode45') subplot(312) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_2') legend('改進(jìn)歐拉法','四階龍格庫塔法','ode45') subplot(313) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_3') legend('改進(jìn)歐拉法','四階龍格庫塔法','ode45')%% 仿真步長(zhǎng) h=0.9 時(shí) % 參數(shù) h = 0.9; % 時(shí)間步長(zhǎng)% 改進(jìn)歐拉法求解 [T1,X1] = ODE_ImprovedEuler( Hfun,t,h,x0 ); % 改進(jìn)歐拉法求解 [T2,X2] = ODE_RK4( Hfun,t,h,x0 ); % Matlab自帶ode45求解 [T3,X3] = ode45( Hfun,t(1):h:t(2),x0 );% 繪圖對(duì)比 figure subplot(311) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_1') legend('改進(jìn)歐拉法','四階龍格庫塔法','ode45') subplot(312) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_2') legend('改進(jìn)歐拉法','四階龍格庫塔法','ode45') subplot(313) plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1)) xlabel('Time(s)') ylabel('x_3') legend('改進(jìn)歐拉法','四階龍格庫塔法','ode45')

不同時(shí)間步長(zhǎng) $h$ 時(shí)的數(shù)值計(jì)算結(jié)果：

步長(zhǎng) $h = 0.6 s$
步長(zhǎng) $h = 0.9 s$

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

4. 聯(lián)系作者

有Matlab/Simulink方面的技術(shù)問題，歡迎發(fā)送郵件至944077462@qq.com討論。
更多Matlab/Simulink原創(chuàng)資料，歡迎關(guān)注微信公眾號(hào)：Matlab Fans

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

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