四阶龙格库塔法(Runge-Kutta)求解常微分方程的 Matlab程序及案例
文章目录1. 算法2. 程序3. 案例4. 联系作者1. 算法上一篇介绍了显式欧拉法、隐式欧拉法、两步欧拉法和改进欧拉法求解常微分方程初值问题；其中显式欧拉法和隐式欧拉法是一阶算法精度，截断误差为O(h2)O\left( {{h^2}} \right)O(h2)；两步欧拉法和改进欧拉法是二阶算法精度，截断误差为O(h3)O\left( {{h^3}} \right)O(h3)；欧拉法的精度有限、需
文章共8,426字 · 阅读需要大约29分钟
1. 算法
上一篇介绍了显式欧拉法、隐式欧拉法、两步欧拉法和改进欧拉法求解常微分方程初值问题；其中显式欧拉法和隐式欧拉法是一阶算法精度，截断误差为 O ( h 2 ) O\left( {{h^2}} \right) O(h2)；两步欧拉法和改进欧拉法是二阶算法精度，截断误差为 O ( h 3 ) O\left( {{h^3}} \right) O(h3)；欧拉法的精度有限、需要求解步长 h h h很小。本篇介绍求解精度更高的四阶龙格库塔法(Runge-Kutta)，其截断误差为 O ( h 5 ) O\left( {{h^5}} \right) O(h5)。
对于一阶微分方程初值问题：
{ y ˙ = f ( t , y ) y ( t 0 ) = y 0 \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. {y˙=f(t,y)y(t0)=y0
式中， t 0 {t_0} t0为初始时间（已知常数）， y 0 {{\bf{y}}_0} y0为初始状态（已知向量）， f ( t , y ) f\left( {t,{\bf{y}}} \right) f(t,y)为关于时间 t {t} t和状态 y {{\bf{y}}} y的函数（已知函数）。
四阶龙格库塔法(Runge-Kutta)求解算法为：
k
1
=
f
(
t
(
k
)
,
y
(
k
)
)
{{k}_{1}}=f\left( t\left( k \right),\mathbf{y}\left( k \right) \right)
k1=f(t(k),y(k))
k
2
=
f
(
t
(
k
)
+
h
2
,
y
(
k
)
+
h
2
k
1
)
{{k}_{2}}=f\left( t\left( k \right)+\frac{h}{2},\mathbf{y}\left( k \right)+\frac{h}{2}{{k}_{1}} \right)
k2=f(t(k)+2h,y(k)+2hk1)
k
3
=
f
(
t
(
k
)
+
h
2
,
y
(
k
)
+
h
2
k
2
)
{{k}_{3}}=f\left( t\left( k \right)+\frac{h}{2},\mathbf{y}\left( k \right)+\frac{h}{2}{{k}_{2}} \right)
k3=f(t(k)+2h,y(k)+2hk2)
k
4
=
f
(
t
(
k
)
+
h
,
y
(
k
)
+
h
k
3
)
{{k}_{4}}=f\left( t\left( k \right)+h,\mathbf{y}\left( k \right)+h{{k}_{3}} \right)
k4=f(t(k)+h,y(k)+hk3)
y
(
k
+
1
)
=
y
(
k
)
+
h
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
\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(k+1)=y(k)+6h(k1+2k2+2k3+k4)
y
(
0
)
=
y
0
\mathbf{y}\left( 0 \right)={{\mathbf{y}}_{0}}
y(0)=y0
上式是关于 y ( k ) {\bf{y}}\left( k \right) y(k)向 y ( k + 1 ) {\bf{y}}\left( k+1 \right) y(k+1)的递推形式，可以根据初始条件按照递推关系依次求出 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(1),y(2),y(3),y(4)⋯,y(N)⋯，此离散序列即为微分方程的数值解。
微分方程的数值解法，本质是使用数值积分来实现
y
˙
{\bf{\dot y}}
y˙向
y
{\bf{y}}
y的转换。四阶龙格库塔法通过对微分的四步分段逼近，在一个求解步长内能够逼近复杂的曲线，因此能够取得较高的计算精度，其截断误差为
O
(
h
5
)
O\left( {{h^5}} \right)
O(h5)。
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}}
x等价于前文中的向量
y
{\bf{y}}
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˙=f(t,x)=⎣⎡x(2)∗x(3)−x(1)∗x(3)−0.51∗x(1)∗x(2)⎦⎤
x
(
0
)
=
[
0
1
1
]
\mathbf{x}\left( 0 \right)=\left[ \begin{matrix} 0 \\ 1 \\ 1 \\ \end{matrix} \right]
x(0)=⎣⎡011⎦⎤
% 四阶龙格库塔算法(Runge-Kutta)测试程序
% By ZFS@wust 2021
% 获取更多Matlab/Simulink原创资料和程序，清关注微信公众号：Matlab Fans
clear
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 h时的数值计算结果：
-
步长 h = 0.6 s h=0.6s h=0.6s
-
步长 h = 0.9 s h=0.9s h=0.9s
从计算结果可以看出，四阶龙格库塔法(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获取百度网盘下载链接。
旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。
更多推荐
57
0- 0
已为社区贡献8条内容
所有评论(0)