四阶龙格库塔方法求解一次常微分方程组

Created2021-06-11Updated2024-05-30

一、写在前面

龙格库塔方法是数值求解常微分非线性方程的有利工具，计算精度较高，通过缩短步进距离和增加阶数可以进一步控制误差范围。工程上较为常用的是四阶龙格库塔算法（R-K4)，在计算收敛的情况下往往可以得到比较好的结果。

二、四阶龙格库塔方法

这里简单介绍一下算法的具体实现过程，不做详细的推导。其求解的问题是形如方程：

y˙=f(y,t),andt∈[t0,t1]inity(t0)=c0\dot{y}=f(y,t), and t\in[t_0,t_1] \\ init y(t_0)=c_0

通过选取一定的步进长度h，来对区间上函数值进行单步迭代求解，最终得到结果。具体计算公式为：

tn+1=tn+hk1=f(yn,tn)k2=f(yn+h2k1,tn+h2)k3=f(yn+h2k2,tn+h2)k4=f(yn+hk3,tn+h)yn+1=yn+h6(k1+2k2+2k3+k4)t_{n+1}=t_n+h\\ k_1=f(y_n,t_n)\\ k_2=f(y_n+\dfrac{h}{2}k_1,t_n+\dfrac{h}{2})\\ k_3=f(y_n+\dfrac{h}{2}k_2,t_n+\dfrac{h}{2})\\ k_4=f(y_n+hk_3,t_n+h)\\ y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)

通过以上公式，选取合适的步进长度h，反复迭代，就可求解出y的数值解。

三、使用四阶龙格库塔方法求解一次常微分方程组

  • 使用的环境matlabR2019a
  • 求解的方程

x˙=y+3z+sin(5t)y˙=x+cos(t)z˙=x+z−3cos(3t)sin(4t)andt∈[0,1]x(0)=y(0)=z(0)=1\dot{x}=y+3z+sin(5t)\\ \dot{y}=x+cos(t)\\ \dot{z}=x+z-3cos(3t)sin(4t)\\ and t\in[0,1]\\ x(0)=y(0)=z(0)=1

  • matlab代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
clear;
clc;
close all;

h=1e-5;      %步进长度
t=0:h:1;   %生成自变量t的向量

%%创建计算结果x，y，z的数组
N=length(t);
x=ones(1,N);
y=ones(1,N);
z=ones(1,N);

%%四阶龙格库塔迭代
for i=2:N
    t_n=t(i-1);
    x_n=x(i-1);
    y_n=y(i-1);
    z_n=z(i-1);
    
    kx1=y_n+3*z_n+sin(5*t_n);
    ky1=x_n+cos(t_n);
    kz1=x_n+z_n-3*cos(3*t_n)*sin(4*t_n);
    
    kx2=(y_n+ky1*h/2)+3*(z_n+kz1*h/2)+sin(5*(t_n+h/2));
    ky2=(x_n+kx1*h/2)+cos(t_n+h/2);
    kz2=(x_n+kx1*h/2)+(z_n+kz1*h/2)-3*cos(3*(t_n+h/2))*sin(4*(t_n+h/2));
    
    kx3=(y_n+ky2*h/2)+3*(z_n+kz2*h/2)+sin(5*(t_n+h/2));
    ky3=(x_n+kx2*h/2)+cos(t_n+h/2);
    kz3=(x_n+kx2*h/2)+(z_n+kz2*h/2)-3*cos(3*(t_n+h/2))*sin(4*(t_n+h/2));
    
    kx4=(y_n+ky3*h)+3*(z_n+kz3*h)+sin(5*(t_n+h));
    ky4=(x_n+kx3*h)+cos(t_n+h);
    kz4=(x_n+kx3*h)+(z_n+kz3*h)-3*cos(3*(t_n+h))*sin(4*(t_n+h));
    
    x(i)=x_n+h/6*(kx1+2*kx2+2*kx3+kx4);
    y(i)=y_n+h/6*(ky1+2*ky2+2*ky3+ky4);
    z(i)=z_n+h/6*(kz1+2*kz2+2*kz3+kz4);  
end
%%画图
figure();
hold on;
plot(t,x,'r');
plot(t,y,'g');
plot(t,z,'b');
legend('x','y','z');
xlabel('t');
title('xyz函数图象');
hold off;

  • 最后得到计算结果图象
    在这里插入图片描述