欧拉方法 - MatLab 中的弹簧振动答案

【问题标题】：Euler's Method - Spring Oscillations in MatLab欧拉方法 - MatLab 中的弹簧振动
【发布时间】：2011-01-27 16:20:27
【问题描述】：

如果我有一个描述质量-弹簧-阻尼系统的一阶微分方程的二阶微分方程，当我不知道第一个方程时，如何使用欧拉法绘制这个方程-阶差分？我想在 MatLab 中执行此操作。这是一个家庭作业问题，所以我没有发布任何代码。我只想简要说明您打算如何做。

二阶导数是 d2(t+1) = (-1/m)*(c*d1+k*y) 其中 c, m, k 是常数，y 最初是 1，d1 是第一个从 0 开始的阶微分，t 是时间。

有什么想法吗？

谢谢:)。

【问题讨论】：

标签： matlab

【解决方案1】：

二阶方程可以转换为一阶微分方程组。

function dy = ex(y)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -c/m*y(2) - k/m*y(1);

从此您可以使用 Matlab 的内置求解器。 ode23s 可以正常工作：

[t,y] = ode23s(@ex, y0, tspan)

【讨论】：

【解决方案2】：

将二阶方程重写为一阶方程组。未知数对应位置和速度。

【讨论】：

【解决方案3】：

一阶微分方程系统可能看起来像这样，而y1 = y。 '用于表示时间导数。

y1' = y2
y2' = -c/m*y2 - k/m*y1

【讨论】：

谢谢。 y2 首先是如何定义的？我想把这些方程放在一个循环中，并在我完成迭代后绘制它们...... :)
@imperial：实现可以在 JudoWill 的帖子中看到。当我们使用二阶系统时，我们需要两个状态变量。 'y2' 的引入是将微分方程重写为一阶微分方程系统的解决方案/模式的一部分。

【解决方案4】：

有一个形式为 (x(t),x'(t)) = exp(A t)*(x(0), x'(0)) 的解析解，其中 A 是一个 2x2 矩阵。如果您不需要使用 matlab 的 ODE 求解器，这是查找系统时间演化的正确方法。

要找到 A，请以这种形式编写您的系统

X = (x, x')
dX = A * X
X = expm(A t) X(0) %注意，要取矩阵的指数，你需要使用 expm，否则你只需单独获取每个元素的指数

所以这里 A = [0 1; -k/m -c/m]

设置k/m = 1，c/m = 0.1，我们可以写

t = linspace(0,20, 1000);
A = [0 1; -1 -0.1];
for j = 1:length(t); x(j,:) = expm([0 1; -1 -0.1]*t(j))*[1;0]; end
plot (t,x)
legend ('position', 'velocity');
title ('underdamped spring starting at y = 1; y'' = 0')

【讨论】：