求解二阶 ODE，Matlab-方程中的加速度需要它自己的值才能包含另一个不同的项

Question

我要在 Matlab 中求解这个二阶 ODE：

(a + f(t))·(dx/dt)·(d²x/dt²)  +  g(t)  +  ((h(t) + i(t)·(d²x/dt² > b·(c-x)))·(dx/dt)  +  j(t))·(dx/dt)²  +  k(t)·(t > d)  = 0

哪里

a、b、c、d 是已知常数
f(t)、g(t)、h(t)、i(t)、j(t)、k(t) 是依赖于 t 的已知函数
x是位置
dx/dt是速度
d²x/dt²是加速度

并注意

的两个条件

(d²x/dt² > b·(c-x))

i(t) 如果 (t > d)
k(t)

所以，这个问题可以在 Matlab 中使用类似的结构来解决，例如：

[T,Y] = ode45(@(t,y) [y(2); 'the expression of the acceleration'], tspan, [x0 v0]);

哪里

T是时间向量，Y是位置向量（第1列为y(1)）和速度向量（第2列为y(2)）。
ode45 是 ODE 求解器，但也可以使用另一个求解器。
tspan,x0,v0 已知。
the expression of the acceleration 表示 d²x/dt² 的表达式，但是 问题来了，因为它在 i(t) 和 'outside' 的条件内同时乘以(a + f(t))·(dx/dt)。所以，加速度不能在matlab中写成d²x/dt² = something

一些可能有帮助的问题：

一旦满足条件(d²x/dt² > b·(c-x))and/or(t > d)，相应的项i(t)and/ork(t)将被介绍到tspan.
对于条件(d²x/dt² > b·(c-x))，项d²x/dt²可以写成速度差，如y(2) - y(2)'，如果y(2)'是前一瞬间的速度除以 tspan 中定义的步进时间。但我不知道如何在求解 ODE

提前谢谢你！

Answer 1

首先，您应该reduce your problem to a first-order differential equation，将dx/dt 替换为速度的动力学变量。这是求解 ODE 时无论如何都必须做的事情，这样您就不需要访问以前的速度值。

至于实现你的条件，只需修改你传递给ode45的函数来解决这个问题。为此，您可以使用 ODE 右侧的 d²x/dt²。请记住，尽管 ODE 求解器不喜欢不连续性，因此您可能希望在满足条件后平滑步骤或使用不同的函数重新启动求解器（归功于 Steve）。

Answer 2

第二个条件项 k(t)*(t>d) 应该很容易实现，所以我将忽略它。

我会将你的等式分成两部分：

F1(t,x,x',x'') := (a+f(t))x'x'' + g(t) + (h(t)x'+j(t))x'' + k(t)(t>d),
F2(t,x,x',x'') := F1(t,x,x',x'') + i(t)x'x'',

其中素数表示时间导数。正如

中所建议的

[...] or just restart the solver with a different function

您可以解决 t \in [t0, t1] =: tspan 的 ODE F1。接下来，您将找到第一个时间 tstar，其中 x''> b(c-x) 以及值 x(tstar) 和 x'(tstar)，并使用 t \in [tstar,t1] 求解 F2 =20=]作为起始条件。

说了这么多，正确的实现应该是使用 events, as suggested in .

Answer 3

所以，我应该使用如下所示的东西：

function [value,isterminal,direction] = ODE_events(t,y,b,c)

value = d²x/dt² - b*(c-y(1));   % detect (d²x/dt² > b·(c-x)), HOW DO I WRITE d²x/dt² HERE?
isterminal = 0;                 % continue integration
direction  = 0;                 % zero can be approached in either direction

然后包含在我的颂歌所在的文件 (.m) 中：

refine = 4; % I do not get exactly how this number would affect the results
options = odeset('Events',@ODE_events,'OutputFcn',@odeplot,'OutputSel',1, 'Refine',refine);

[T,Y] = ode45(@(t,y) [y(2); ( 1/(a + f(t))*(y(2)))*( - g(t) - ((h(t) + i(t))·(y(2))  -  j(t)·(y(2))² - k(t)*(t > d))  ], tspan, [x0 v0], options);

如何处理i(t)？因为必须以某种方式包含 i(t)*(d²x/dt² > b*(c-y(1))))*(y(2))。

再次感谢您

求解二阶 ODE，Matlab-方程中的加速度需要它自己的值才能包含另一个不同的项

Solving 2nd order ODE, Matlab- the acceleration in the equation needs its own value in order to include another different term

matlab

position

acceleration

ode