Решая ОДУ 2-го порядка, Matlab - ускорение в уравнении нуждается в собственном значении, чтобы включить еще один другой член

У меня есть ODE 2-го порядка для решения в 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

куда

  • _2 _, _ 3 _, _ 4 _, _ 5_ - известные константы
  • _6 _, _ 7 _, _ 8 _, _ 9 _, _ 10 _, _ 11_ - известные функции, зависящие от t
  • x - это позиция
  • dx/dt - скорость
  • d²x/dt² - это ускорение

и обратите внимание на два условия, которые

  • i(t) вводится в уравнение, если (d²x/dt² > b·(c-x))
  • 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, но можно использовать и другой.
  • _26 _, _ 27 _, _ 28_ известны.
  • 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

Заранее благодарю !


matlab position ode acceleration
person Sergio    schedule 24.09.2017    source источник
comment
Вы можете попробовать неявный решатель, например ode15i   -  person Steve    schedule 25.09.2017
comment
@ Стив: И что это решит?   -  person Wrzlprmft    schedule 25.09.2017
comment
@Wrzlprmft Это решает тот факт, что это неявное ODE, где x_tt := d^2x/dt^2 является эффективно неявным аргументом как H(x_tt - b(c*x)), где H(t) - это ступенчатая функция Хевисайда H(t) = 1 if t>=0, 0 if t<0. Существует ли решение или имеет ли он смысл - это еще один вопрос, на который стоит обратить внимание.   -  person Steve    schedule 25.09.2017
comment
На самом деле то, что я сказал, не совсем верно, из-за вашего условия , как только условие [...] будет выполнено, соответствующие [условия] будут введены до конца определенного времени   -  person Steve    schedule 25.09.2017

Ответы (3)


arrow_upward
2
arrow_downward

Прежде всего, вам следует свести проблему к дифференциальному уравнению первого порядка, выполнив заменяя dx / dt динамической переменной вместо скорости. Это то, что вам все равно нужно делать для решения ODE, и таким образом вам не нужно обращаться к предыдущим значениям скорости.

Что касается реализации ваших условий, просто измените функцию, которую вы передаете ode45, чтобы учесть это. Для этой цели вы можете использовать то, что d²x / dt² находится в правой части вашего ODE. Однако имейте в виду, что решатели ODE не любят разрывов, поэтому вы можете сгладить шаг или просто перезапустить решатель с другой функцией, как только условие будет выполнено (кредит на Стив).

person Wrzlprmft    schedule 25.09.2017
comment
Я полностью согласен с первым абзацем, но относительно второго ode45 полагаю, что вы можете выразить уравнение в терминах производной высшего порядка, что невозможно из-за i(t)*(x_tt>b(c-x))*x_t*x_tt члена - person Steve; 25.09.2017
comment
В любом случае вы не можете явно оценить это состояние, поскольку оно остается после включения, и вам в любом случае придется полагаться на отсутствие парадокса. (Также см. Мою правку.) - person Wrzlprmft; 25.09.2017
comment
+1 Хорошо, теперь это имеет для меня больше смысла. Я начал писать перезапуск решателя с другой функциональной частью, поэтому я закончу и опубликую это, но этот ответ, похоже, покрыл это. - person Steve; 25.09.2017
comment
Вы должны использовать events, чтобы найти точки для перезапуска интеграции, используя 0-1-константы, которые изменяются во время события, вместо того, чтобы иметь локальные условия внутри функции ODE. Таким образом интеграция будет плавной и может надежно найти конечные точки каждого сегмента. - person Lutz Lehmann; 25.09.2017
comment
@JRM: прочтите документацию и попробуйте реализовать это самостоятельно (задайте дополнительный вопрос, если у вас возникнут проблемы). Я не говорю на Matlab. - person Wrzlprmft; 26.09.2017
comment
См. Новый ответ @Wrzlprmft - person Sergio; 30.09.2017
comment
@JRM: Нового ответа нет. Если вы хотите ответить на свой вопрос самостоятельно, это нормально, но, пожалуйста, опубликуйте фактический ответ и не редактируйте свой вопрос. - person Wrzlprmft; 30.09.2017

arrow_upward
1
arrow_downward

Второй условный термин 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'',

где штрихом обозначены производные по времени. Как предлагается в этом другом ответе

[...] или просто перезапустите решатель с другой функцией

вы можете решить 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) в качестве начальных условий.

Сказав все это, для правильной реализации этого следует использовать < em> события, как предлагается в комментарий LutzL.

person Steve    schedule 25.09.2017

arrow_upward
0
arrow_downward

Итак, я должен использовать что-то вроде этого:

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)) нужно как-то включить.

Еще раз спасибо

person Sergio    schedule 30.09.2017