如何使用 ode23 求解这个微分方程？答案

【问题标题】：How to solve this differential equation using ode23?如何使用 ode23 求解这个微分方程？
【发布时间】：2021-02-01 15:54:32
【问题描述】：

我这里有一个齿轮系统，我需要使用它来模拟 MATLAB 和 Simulink。

为此，首先我需要得到这个的微分方程旋转（齿轮）系统。我认为我设法做到了，就像关注

iT1 = J * theta2''+D * theta2'+K * theta2,

其中J = i^2 * Jn1 + Jn2 + J2（另外，i 是齿轮比；i = theta1/theta2 = N2/N1 = Tn2 / Tn1 = T2/T1）。

现在我需要求解这个二阶微分方程 MATLAB 中的 ode23。问题是，我不知道任何初始条件，我也不知道如何从这个特定系统中获取它们，例如y(0) = 0, y'(0) = 1 如果没有初始条件，还是我要建立初始条件...？如果任何人都需要它我的问题的参数如下：

N2 = 90, N1 = 36, 
Jn1 = 0.5 Nms^2/rad, Jn2 = 0.8 Nms^2/rad, 
J2 = 2 Nms^2/rad, D = 8 Nms/rad, K = 5 Nm/rad

y=theta2 和 x=T1

    y''+1.3502y'+0.8439y = 0.4219x

任何帮助将不胜感激！

【问题讨论】：

您能多描述一下机械系统吗？一侧有齿轮组件，另一侧有阻尼弹簧，它们是如何连接的？您是否正在寻找像极限环或其他一些渐近线这样的长期行为，其中任何初始条件的影响都变得可以忽略不计？ // 这里的问题可能是题外话，scicomp.SE 或physics.SE 可能更合适，检查它们的范围描述。理想情况下，您必须在此处提供一些您遇到问题的代码。
我把这个特殊的系统想象成一个齿轮箱，第二个轴上有一个飞轮。我想要做的是模拟这个特定系统的输出，用于步进输入和正弦波输入。我的输出应该是第二个轴的位置，以弧度为单位。

标签： matlab system simulink ode

【解决方案1】：

因为我不知道您对微分方程的教育，首先：正如我的教授一直所说的那样，初始条件或边界条件（在这种情况下为初始条件）构成了解决方案，因为任何 ode 都通过积分求解，没有这些条件，实际上无法求解 ode。

但是看起来你似乎有一个二阶颂歌，所以你可以确定一个描述系统本身的解决方案。

我不知道如何在这里写乳胶，但一个简单的例子：

给定最简单的颂歌

dy/dx = y

这可以通过变量的拆分来解决

1/y * dy = dx

集成

ln(y) = x + C_0 -> y = exp(x) + C 其中 C 为 exp(c_0)

您的 ode 应该可以做到这一点。查看 matlab 中的符号数学工具箱或阅读有关 odes 的信息（我没有时间，但在我看来，您的 ode 应该可以解析求解 - 搜索：解析求解二阶微分方程）。

关于 C 的一些话，实际上这是你的边界或初始条件结束的地方。对于上面的例子，这可能是人口的增长，那么 C 将是您的初始人口规模。如果没有初始条件，您只能说系统将如何演变，而不能说在特定时间点它将处于哪种状态。

另一个词：其实你对我的颂歌似乎是众所周知的谐波振荡器。您会为此找到很多东西。

也许您想确定系统的特征值以确定对于任何频率是否稳定，那么您的方程必须是 x'' + dx' + kx = 0 的形式（必须为零） - 这将最终解决我的示例中的特征方程（lambda 中的二阶多项式）

如果您最终获得初始值，则要通过 ode23 求解，您必须将其设为微分方程组（谷歌搜索或查看数学书），然后求解该系统。在 MATLAB 中编写一个函数，该函数返回微分方程组的梯度，例如

function grads = calc_grads(immediate_vals)

% calculate all the n grads - here two - determined by the ode
grads(1) = ...
grads(2) = ...
end

此函数必须可由 ode23 调用。

【讨论】：

【解决方案2】：

代码应该是这样的

% parameters
N2  = 90;
N1  = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2  = 2;
D   = 8;
K   = 5;
J   = (N2/N1)^2 * Jn1 + Jn2 + J2;

% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);

% initial state: (position, velocity) [rad; rad/s]
x0 = [0.1; 0];

% define the time span
t = linspace(0, 15, 10000)';

% define the input step
T1 = zeros(length(t), 1);
T1(t>=1) = 1;

% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);

% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);

% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
       'location', 'southeast');
hold off;

我现在想知道的是如何模拟正弦波函数的输出？有什么帮助吗？

【讨论】：

【解决方案3】：

根据文档，“MATLAB ODE 求解器仅求解一阶方程。您必须将高阶 ODE 重写为一阶方程的等效系统。”[1] 也就是说，要求解第二个方程阶方程如：

iT1 = J*theta2'' + D*theta2' + K*theta2

你应该把它改写成：[1,2]

Y = [  Y_1  ] = [ theta2  ]
    [  Y_2  ]   [ theta2' ]

d/dt Y = [             Y_2              ]
         [ -1/J * (D*Y_2 + K*Y_1 - iT1) ]

其中t 是自变量。 d/dt Y_1 = Y_2 的定义可能看起来是重言式，但有必要统一其他不相关的（就机器所知）变量 Y_1 和 Y_2。

关于初始值，是的，它们是必要的；你必须知道向量Y的初始值。

请记住，ode23() 是一种数值求解器（MATLAB 提供的众多解决方案之一），可用于许多高度复杂的方程，尤其是那些无法解析求解的方程。如果您想要不需要预先固定初始条件的分析解决方案，ode23() 不是适合这项工作的工具；如上一个答案中突出显示的符号数学工具箱在这种情况下可能更有用。然而，迟早你可能最终需要定义初始条件。详见[3]。

[1] https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html（参见“高阶 ODE”部分）

[2]https://en.wikipedia.org/wiki/Ordinary_differential_equation#Reduction_of_order

[3]https://en.wikipedia.org/wiki/Initial_value_problem

【讨论】：