# How to solve a system of second order ODE from time t = T to t = 0

I have a system of second-order ODEs
$$\mathbf{M\ddot{x} + C\dot{x} + Kx = f}$$ I want to know some good numerical methods to solve this system of the equation given the initial conditions at time $$t = T$$ that is I am given $$x(T)$$ and $$\dot{x}(T)$$ and I want to solve for $$x$$ between $$t = [0, T)$$
I have tried to solve this with the Newmark method which is stable with a positive time step but with a negative time step, Newmark fails miserably.

• Does $M$, $C$, $K$ or $f$ vary with time? – fibonatic Apr 11 at 21:52
• Yes K and f varies with time – Random Techie Apr 12 at 14:52

This has a direct analytical solution that in probably most softwares can be solved in one line...

For example mathematica:

DSolve[{m x''[t] + d x'[t] + c x[t] == f, x'[T] == a, x[T] == b}, x[t], t]


Or, if one wanted...by hand using pretty traditional ODE methods.

• Does this give you anything useful when the model has more than a single degree of freedom? It's not very obvious what the expression even mean, if m c and d are square matrices. – alephzero Apr 11 at 23:47
• Ahh, i apparently didn‘t pay attention, i would need to see the system, this solution is only for one dof...a coupled system doesn‘t necessarily have an analytical solution no, though depending on its order it may still. – morbo Apr 12 at 7:09

Your initial second order ODE can also be written as the following system of first order ODE

$$\begin{bmatrix} \dot{x} \\ \ddot{x} \end{bmatrix} = \begin{bmatrix} 0 & I \\ -M^{-1} K & -M^{-1} C \end{bmatrix} \begin{bmatrix} x \\ \dot{x} \end{bmatrix} + \begin{bmatrix} 0 \\ M^{-1} f \end{bmatrix}, \tag{1}$$

which can be solved by many known solvers. The terminal condition can also be used as an initial condition when also negating the right-hand-side of $$(1)$$. Namely, given that

$$\frac{dx}{dt} = f(t,x), \quad x(T)=x_T. \tag{2}$$

Now by using that $$\tau=T-t$$ and $$z(\tau)=x(T-\tau)$$ then $$(2)$$ would be equivalent to

$$\frac{dz}{d\tau} = \frac{dx}{dt} \frac{dt}{d\tau} = -f(T-\tau,z), \quad z(0)=x_T. \tag{3}$$

However, it might be possible that the negated dynamics is unstable, which potentially causes the solution backwards in time to diverge. This would mean that also the analytical solution might also diverge.

• Thank you for the idea. The conversion of the second-order system to the first order is pretty clear. – Random Techie Apr 12 at 15:05

The problem with stability is caused by the equation itself, not by your attempts at solving it.

If you integrate backwards, the "damping terms" are adding energy to the system as the integration proceeds, not dissipating it. Any numerical rounding errors will therefore grow with time.

Even if the numerical errors are negligible, the initial conditions are very poorly determined by the final state of the system.

• Thank you for your insight. So if we consider zero damping and my forward integration is stable then backward integration should also be stable, right? – Random Techie Apr 12 at 15:09