用四阶RK算法编程计算求解简单的振动微分方程并画出曲线

使用四阶龙格-库塔（Runge-Kutta）方法来求解一个简单的振动微分方程，并绘制其解的曲线。假设我们要解的振动微分方程是：

\(\ddot{x} + \omega^2 x = 0\)

这是一个简谐振动方程，其中 (\omega) 是角频率。为了使用四阶龙格-库塔方法，我们需要将其转换为一阶微分方程组。设：

\(x_1 = x\)
\(x_2 = \dot{x}\)

则原方程可以写成：

\(\dot{x_1} = x_2\)
\(\dot{x_2} = -\omega^2 x_1\)

我们用四阶龙格-库塔方法求解这个方程组，并用MATLAB绘制解的曲线。

代码

% 参数设置
omega = 1; % 角频率
t0 = 0; % 初始时间
tf = 10; % 终止时间
h = 0.01; % 时间步长
N = (tf - t0) / h; % 时间步数

% 初始条件
x1_0 = 1; % 初始位移
x2_0 = 0; % 初始速度

% 初始化解数组
t = t0:h:tf;
x1 = zeros(1, N+1);
x2 = zeros(1, N+1);
x1(1) = x1_0;
x2(1) = x2_0;

% 四阶龙格-库塔方法
for i = 1:N
    k1_1 = x2(i);
    k1_2 = -omega^2 * x1(i);
    
    k2_1 = x2(i) + h * k1_2 / 2;
    k2_2 = -omega^2 * (x1(i) + h * k1_1 / 2);
    
    k3_1 = x2(i) + h * k2_2 / 2;
    k3_2 = -omega^2 * (x1(i) + h * k2_1 / 2);
    
    k4_1 = x2(i) + h * k3_2;
    k4_2 = -omega^2 * (x1(i) + h * k3_1);
    
    x1(i+1) = x1(i) + (h/6) * (k1_1 + 2*k2_1 + 2*k3_1 + k4_1);
    x2(i+1) = x2(i) + (h/6) * (k1_2 + 2*k2_2 + 2*k3_2 + k4_2);
end

% 绘制结果
figure;
plot(t, x1, 'b', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Displacement (x)');
title('Displacement vs. Time');
grid on;

figure;
plot(t, x2, 'r', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Velocity (\dot{x})');
title('Velocity vs. Time');
grid on;

说明

1. 参数设置：

   • omega：角频率。
   • t0 和 tf：时间范围。
   • h：时间步长。
   • N：总的时间步数。

2. 初始条件：

   • x1_0：初始位移。
   • x2_0：初始速度。

3. 四阶龙格-库塔方法：

   • 对于每个时间步，计算四个中间值 \(k1, k2, k3, k4\)，然后根据这些中间值更新 \(x_1\)和 \(x_2\)。

4. 绘图：

   • 绘制位移 \(x_1\) 和速度 \(x_2\) 随时间的变化曲线。

参考仿真代码 用四阶RK算法编程计算求解简单的振动微分方程并画出曲线 www.youwenfan.com/contentcnd/99547.html

运行结果

运行上述代码后，你会得到两个图：

1. 位移随时间的变化：显示简谐振动的位移曲线。
2. 速度随时间的变化：显示简谐振动的速度曲线。