解决二阶微分方程，Matlab-方程中的加速需求，以包括其他不同期限

问题描述：

我有这样的二阶微分方程在Matlab解决其自身的价值：解决二阶微分方程，Matlab-方程中的加速需求，以包括其他不同期限

(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)，依赖t
xj(t)，k(t)是已知的功能是位置
dx/dt是速度
d²x/dt²是加速度

并注意两个条件

如果引入
k(t)是在方程中引入如果(t > d)

所以，这个问题可以用类似的结构在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)，并在同一时间“外部”乘以(a + f(t))·(dx/dt)的条件内。因此，加速度不能在MATLAB写为d²x/dt² = something

的一些问题，可以帮助：

一次(d²x/dt² > b·(c-x))和/或(t > d)满足的条件下，各期限i(t)和/或k(t)会被引入到tspan的确定时间的末尾。
为条件(d²x/dt² > b·(c-x))，术语d²x/dt²可以写成速度的差，像y(2) - y(2)'，如果y(2)'是前瞬间的速度，通过在tspan限定的步进时间划分。但我不知道如何在访问速度的前值的ODE

谢谢你在先进的解决！

你可以尝试隐式求解器，如['ode15i'（https://uk.mathworks.com/help/matlab/ref/ode15i.html） – Steve

@Steve：和那会解决什么？ – Wrzlprmft

@Wrzlprmft这解决了这是一个隐式ODE，其中'x_tt：= d^2x/dt^2 *有效*隐式参数为'H（x_tt-b（c * x））'，其中'如果t> = 0，则H（t）是Heaviside阶跃函数，H（t）= 1，如果t Steve

答

首先，你应该reduce your problem to a first-order differential equation，用动力学变量代替速度的dx/dt。这是你必须要做的事情来解决ODE，这种方式你不需要访问以前的速度值。

至于实现您的条件，只需修改您传递给ode45的功能来解决这个问题。为此，您可以使用d²x/dt²位于ODE的右侧。请记住，虽然ODE求解器不喜欢不连续性，所以一旦满足条件（信用额度为Steve），您可能希望平滑该步骤或仅使用不同函数重新启动求解器。

我完全同意第一段，但关于第二段，'ode45'依赖于你可以用最高阶导数表达方程，这是由于'i（t）*（x_tt> b（cx））* x_t * x_tt'术语 – Steve

无论如何，您无法明确评估这种情况，因为它一旦打开就无法继续使用，您必须依赖于悖论。（另见我的编辑。） – Wrzlprmft

+1好吧，现在对我更有意义。我已经开始编写“用不同的函数重新启动解算器”部分，所以我会完成并发布它，但是这个答案似乎已经涵盖了它。 – Steve

答

第二个条件术语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'',

，其中黄金是时间轴衍生物。作为this other answer

建议[...]或只需重新启动求解器具有不同的功能

你能解决ODE F1为t \in [t0, t1] =: tspan。接下来，您会首次发现tstar，其中x''> b(c-x)和值x(tstar)和x'(tstar)，并且解决F2对于t \in [tstar,t1]和x(tstar), x'(tstar)作为起始条件。

说了这么多，正确执行此操作应该使用events，如LutzL's comment中所建议的。

答

所以，我应该用一些看起来是这样的：

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))必须以某种方式包含。

再次感谢您

解决二阶微分方程，Matlab-方程中的加速需求，以包括其他不同期限

相关推荐