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

Question

I have this 2nd order ODE to solve in 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

where

• a,b,c,d are known constants
• f(t),g(t),h(t),i(t),j(t),k(t) are known functions dependent on t
• x is the position
• dx/dt is the velocity
• d²x/dt² is the acceleration

and notice the two conditions that

• i(t) is introduced in the equation if (d²x/dt² > b·(c-x))
• k(t) is introduced in the equation if (t > d)

So, the problem could be solved with a similar structure in Matlab as this example:

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

where

• T is the time vector, Y is the vector of position (column 1 as y(1)) and velocity (column 2 as y(2)).
• ode45 is the ODE solver, but another one could be used.
• tspan,x0,v0 are known.
• the expression of the acceleration means an expression for d²x/dt², but here comes the problem, since it is inside the condition for i(t) and 'outside' at the same time multiplying (a + f(t))·(dx/dt). So, the acceleration cannot be written in matlab as d²x/dt² = something

Some issues that could help:

• once the condition (d²x/dt² > b·(c-x)) and/or (t > d) is satisfied, the respective term i(t) and/or k(t) will be introduced until the end of the determined time in tspan.

• for the condition (d²x/dt² > b·(c-x)), the term d²x/dt² could be written as the difference of velocities, like y(2) - y(2)', if y(2)' is the velocity of the previous instant, divided by the step-time defined in tspan. But I do not know how to access the previous value of the velocity during the solving of the ODE

Thank you in advanced !

Answer 1

First of all, you should reduce your problem to a first-order differential equation, by substituting dx/dt with a dynamical variable for the velocity. This is something you have to do anyway for solving the ODE and this way you do not need to access the previous values of the velocity.

As for realising your conditions, just modify the function you pass to ode45 to account for this. For this purpose you can use that d²x/dt² is in the right-hand side of your ODE. Keep in mind though that ODE solvers do not like discontinuities, so you may want to smoothen the step or just restart the solver with a different function, once the condition is met (credit to Steve).

Answer 2

The second conditional term k(t)*(t>d) should be simple enough to implement, so I'll pass over that.

I would split up your equation into two part:

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'',

where prime denote time derivatives. As suggested in this other answer

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

you could solve the ODE F1 for t \in [t0, t1] =: tspan. Next, you'd find the first time tstar where x''> b(c-x) and the values x(tstar) and x'(tstar), and solve F2 for t \in [tstar,t1] with x(tstar), x'(tstar) as starting conditions.

Having said all that, the proper implementation of this should be using events, as suggested in LutzL's comment.

Answer 3

So, I should use something that looks like this:

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

and then include in the file (.m), where my ode is, this:

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);

How do I handle i(t)? because i(t)*(d²x/dt² > b*(c-y(1))))*(y(2)) has to be included somehow.

Thank you again