解決二階微分方程，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

謝謝你在先進的解決！

來源

2017-09-24 Sergio

你可以嘗試隱式求解器，如['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 <0，則爲0。解決方案是存在還是有意義是另一個需要注意的問題。 – Steve

首先，你應該reduce your problem to a first-order differential equation，用動力學變量代替速度的dx/dt。這是你必須要做的事情來解決ODE，這種方式你不需要訪問以前的速度值。

至於實現您的條件，只需修改您傳遞給ode45的功能來解決這個問題。爲此，您可以使用d²x/dt²位於ODE的右側。請記住，雖然ODE求解器不喜歡不連續性，所以一旦滿足條件（信用額度爲Steve），您可能希望平滑該步驟或僅使用不同函數重新啓動求解器。

來源

2017-09-25 08:58:10 Wrzlprmft

我完全同意第一段，但關於第二段，'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中所建議的。

來源

2017-09-25 09:37:18 Steve

所以，我應該用一些看起來是這樣的：

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

再次感謝您

來源

2017-09-30 09:06:23 Sergio

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

回答

相關問題