如何用Octave求解二阶微分方程组并绘制h-时间曲线？

完全可行。你已经具备了数值求解二阶微分方程组所需的全部必要条件：包括h、θ的初始值，它们的一阶导数初始值，以及所有固定参数（μ、r、g、L）。只要将二阶方程组转化为一阶方程组，就能用Octave的ODE求解器处理，最终绘制h随时间变化的曲线。

1. 二阶方程组转一阶方程组 #

二阶微分方程的标准处理方式是引入新变量拆分阶数。假设你的二阶方程组形式为：

$\ddot{h} = f(h, \dot{h}, \theta, \dot{\theta}, \mu, r, g, L, t)$
$\ddot{\theta} = g(h, \dot{h}, \theta, \dot{\theta}, \mu, r, g, L, t)$

引入状态向量 $y = [h, \dot{h}, \theta, \dot{\theta}]^T$，转化为以下一阶方程组：

• $\dot{y}_1 = y_2$ （对应 $\dot{h} = y_2$）
• $\dot{y}_2 = f(y_1, y_2, y_3, y_4, \mu, r, g, L, t)$ （对应 $\ddot{h}$ 的表达式）
• $\dot{y}_3 = y_4$ （对应 $\dot{\theta} = y_4$）
• $\dot{y}_4 = g(y_1, y_2, y_3, y_4, \mu, r, g, L, t)$ （对应 $\ddot{\theta}$ 的表达式）

2. 编写Octave的ODE函数 #

创建一个名为system_ode.m的文件，填入以下内容（替换为你实际的方程组表达式）：

function dydt = system_ode(t, y, miu, r, g, L)
    % 拆分状态变量
    h = y(1);
    dh_dt = y(2);
    theta = y(3);
    dtheta_dt = y(4);
    
    % 替换为你的二阶导数公式
    d2h_dt2 = ...; % 这里写ddot(h)的具体计算式
    d2theta_dt2 = ...; % 这里写ddot(theta)的具体计算式
    
    % 组装一阶方程组的导数输出
    dydt = [dh_dt;
            d2h_dt2;
            dtheta_dt;
            d2theta_dt2];
end

3. 设置初始条件并求解 #

在Octave命令行或脚本中执行：

% 定义参数（替换为你的实际数值）
miu = 0.1;
r = 0.5;
g = 9.81;
L = 2.0;

% 初始条件：[h初始值, h的一阶导初始值, theta初始值, theta的一阶导初始值]
y0 = [1.0; 0.0; pi/4; 0.0]; % 替换为你的实际初始值

% 求解时间范围，示例为t从0到10
tspan = [0, 10];

% 调用ode45求解，传递额外参数
[t, y] = ode45(@(t,y) system_ode(t,y,miu,r,g,L), tspan, y0);

4. 绘制h-t曲线 #

plot(t, y(:,1), 'LineWidth', 1.5);
xlabel('时间 t');
ylabel('h(t)');
title('h随时间变化曲线');
grid on;

• 检查方程组表达式：确保变量、参数拼写正确，乘法运算符*不能省略（Octave中ab是变量名，a*b才是乘法）
• 初始条件维度：y0必须是4维列向量，与状态向量维度匹配
• ODE函数输出：dydt必须是列向量，和输入状态向量的维度一致
• 参数传递：通过匿名函数@(t,y) system_ode(t,y,miu,r,g,L)确保参数正确传入ODE函数

内容的提问来源于stack exchange，提问作者Zhi Ning