有时候不方便用MATLAB中自带的求解器，写一个定步长的四阶龙格库塔求解器。

四阶龙格库塔function

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
function [t, y] = rungeKutta4(odefun, tspan, y0, h)
    % 输入参数:
    % odefun - 常微分方程的函数句柄 (dy/dt = f(t, y))
    % tspan  - 时间区间 [t0, tf]
    % y0     - 初始条件
    % h      - 时间步长
  
    % 输出参数:
    % t      - 时间点数组
    % y      - 解的数组

    % 初始化
    t0 = tspan(1);
    tf = tspan(2);
    t = t0:h:tf;  % 生成时间点
    n = length(t);
    y = zeros(length(y0), n);
    y(:,1) = y0;
  
    % 四阶龙格库塔算法
    for i = 1:n-1
        k1 = h * odefun(t(i), y(:,i));
        k2 = h * odefun(t(i) + 0.5*h, y(:,i) + 0.5*k1);
        k3 = h * odefun(t(i) + 0.5*h, y(:,i) + 0.5*k2);
        k4 = h * odefun(t(i) + h, y(:,i) + k3);
  
        y(:,i+1) = y(:,i) + (k1 + 2*k2 + 2*k3 + k4) / 6;
    end
end

二阶弹簧阻尼系统

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 function dx = damped_system(t, x)
    % 输入参数：
    %   t - 时间
    %   x - 状态向量，x(1) 为位移 x，x(2) 为速度 v = dx/dt
    %   m - 质量
    %   c - 阻尼系数
    %   k - 弹簧常数
    % 输出：
    %   dx - 状态向量的导数，dx(1) 为速度，dx(2) 为加速度
    m=2;
    c=0.9;
    k=0.5;
    % 初始化导数向量
    dx = zeros(2, 1);
  
    % 一阶导数为速度
    dx(1) = x(2);
  
    % 二阶导数为加速度，根据二阶阻尼系统方程
    dx(2) = -(c/m) * x(2) - (k/m) * x(1);
end

使用RK4求解二阶弹簧阻尼

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
% 定义微分方程
odefun = @damped_system;

% 定义时间区间、初始条件和时间步长
tspan = [0, 20];  % 从0到2秒的时间区间
y0 = [0.1; 0];     % 初始条件，x1(0) = 1, x2(0) = 1
h = 0.1;         % 时间步长

% 使用四阶龙格库塔求解器
[t, y] = rungeKutta4(odefun, tspan, y0, h);

% 绘制结果
figure;
subplot(2, 1, 1);
plot(t, y(1, :), '-o');
xlabel('Time t');
ylabel('x1(t)');
title('Solution of x1(t) using Runge-Kutta 4th order method');
grid on;

subplot(2, 1, 2);
plot(t, y(2, :), '-o');
xlabel('Time t');
ylabel('x2(t)');
title('Solution of x2(t) using Runge-Kutta 4th order method');
grid on;

运行结果：

rk4