# Decouple 3 simultaneous 2nd order ODEs (of 3 dependent variables)

I have three simultaneous ODEs with coupled dependent variables. I was wondering whether there's a way for Mathematica to decouple these equations (into three ODEs of a single dependent variable each), rather than attempt to solve for the variables (which it currently fails to do).

Specifically, the (messy) equations are $$\begin{cases} k x_k + k x_A(t) + m \ddot{x}_A(t) = F_D(t) + k x_B(t) \\ 2 k x_B(t) + m \ddot{x}_B(t) = k(x_A(t) + x_C(t)) \\ m \ddot{x}_C(t) = k(x_k + x_B(t) - x_C(t)) \end{cases}$$ where everything that isn't $F_D(t), x_A(t), x_B(t), x_C(t)$ or their derivatives is a constant ($m, k, x_k$). Accompanying is a bunch of initial conditions: $$x_A(0) = 0, \;x_B(0) = x_k, \;x_C(0) = 2 x_k$$ $$\dot{x}_A(0) = \dot{x}_B(0) = \dot{x}_C(0) = 0$$ $$\ddot{x}_A(0) = \ddot{x}_B(0) = \ddot{x}_C(0) = 0$$

From this, I'd like Mathematica to return an equation of the form $$\ddot{x}_C(t) + \alpha \dot{x}_C(t) + \beta x_C(t) = f(F_D(t))$$ (I'm only interested in one of the 3 ODEs of this form).

I should note I tried unsuccessfully to use DSolve, which ran for a very long time then reported For some branches of the general solution, unable to compute the limit at the given points. Some of the solutions may be lost. and output {}.

My input was (letting $F_D(t) = F_0 \cos(t \omega)$)

DSolve[{k xk + k xa[t] + m xa''[t] ==
F0 Cos[t \[Omega]] + k xb[t],
2 k xb[t] + m xb''[t] == k (xa[t] + xc[t]),
m xc'')[t] == k (xk + xb[t] - xc[t]), xa'[0] == 0,
xb'[0] == 0, xc'[0] == 0, xa''[0] ==  0, xb''[0] ==  0,
xc''[0] ==  0, xa[0] == 0, xb[0] ==  xk, xc[0] == 2 xk}, {xa[t],
xb[t], xc[t]}, t].


So, is there anyway to have Mathematica decouple these equations into isolated variables without attempting to solve for the dependent variables? I've heard something vaguely related about using matrices.

Thanks!

• Why impose initial values for the second derivatives if your ODE is second-order? – J. M.'s discontentment Aug 30 '15 at 12:04
• They're only superfluous since I already have 6 initial conditions; they're a physical result of the first presented conditions. I thought it might make Mathematica's job a little easier. – Anti Earth Aug 30 '15 at 12:06
• The system is linear with constant coefficients and with a harmonic inhomogeneity. You can write it as a Matrix ODE for the vector x of functions with: x'' = M.x + f, and proceed as if it were a scalar ODE using MMA's function MatrixExp. Sorry, only GedankenMMA available to me on sunday. – Dr. Wolfgang Hintze Aug 30 '15 at 12:18
• I forgot to mention: decoupling the equations should be doable using EigenSystem. – Dr. Wolfgang Hintze Aug 30 '15 at 12:24
• Can you point me to anymore information? I'm not familiar with solving systems of in-homogeneous equations with matrices. – Anti Earth Aug 31 '15 at 12:10

This is a nice exercise.

First of all, we write down the system of equations in vector form as follows (letting m->1, k->1 for simplicity, and without loss of generality)

vx''[t] == ma.vx[t] + vf[t];


where

vx[t_] = {xA[t], xB[t], xC[t]};

ma = {{-1, 1, 0}, {1, -2, 1}, {0, 1, -1}};


and

vf[t_] = {FD[t] - xk, 0, xk};


Check that the vector equation is reproduced:

eq = D[vx[t], {t, 2}] == ma.vx[t] + vf[t] // Thread

(*
Out[5]= {(xA^\[Prime]\[Prime])[t] == -xk + FD[t] - xA[t] + xB[t], (
xB^\[Prime]\[Prime])[t] ==
xA[t] - 2 xB[t] + xC[t], (xC^\[Prime]\[Prime])[t] == xk + xB[t]-xC[t]}
*)


In order to decouple the equations we look for the Eigensystem of the matrix ma

es = Eigensystem[ma]

(*
Out[6]= {{-3, -1, 0}, {{1, -2, 1}, {-1, 0, 1}, {1, 1, 1}}}
*)


The first component consists of the eigenvalues, the second is the matrix of the corresponding eigenvectors.

Let the letter matrix be

u = es[[2]]

(*
Out[7]= {{1, -2, 1}, {-1, 0, 1}, {1, 1, 1}}
*)


This matrix product

md = u.ma.Inverse[u]

(*
Out[8]= {{-3, 0, 0}, {0, -1, 0}, {0, 0, 0}}
*)


gives a diagonal matrix with the eigenvalues in the main diagonal.

Now matrix - multipy eq from the left by u. Because u is independent of time we can write

D[u.vx[t], {t, 2}] == u.ma.Inverse[u].(u.vx[t]) + u.vf[t];


Letting

vg[t_] = u.vf[t];


we obtain a vector equation for the vector u.vx[t] which we call vz[t].

vz[t_] = {vz1[t], vz2[t], vz3[t]}

(*
Out[18]= {vz1[t], vz2[t], vz3[t]}
*)

eq1 = D[vz[t], {t, 2}] == md.vz[t] + vg[t] // Thread;
Column[%]


$$\left( \begin{array}{c} \text{vz1}''(t)=\text{FD}(t)-3 \text{vz1}(t) \\ \text{vz2}''(t)=2 \text{xk}-\text{FD}(t)-\text{vz2}(t) \\ \text{vz3}''(t)=\text{FD}(t) \\ \end{array} \right)$$

Hence we have three equations, each one only for one variable (component of vector vz)

The general solution of the homogeneous equation (xk->0, FD->0) is

solh = DSolve[eq1, vz[t], t] /. {xk -> 0, FD[_] -> 0}

(*
Out[35]= {{vz1[t] -> C[1] Cos[Sqrt[3] t] + C[2] Sin[Sqrt[3] t],
vz2[t] -> C[3] Cos[t] + C[4] Sin[t], vz3[t] -> C[5] + t C[6]}}
*)


As it should it contains 6 integration constants to be determined by the initial conditions.

Before doing this, we need to transform back from the vector vz to the original vector vx, called vxs

This is done by inverting

vxs[t_] = Inverse[u].vz[t]

(*
Out[40]= {vz1[t]/6 - vz2[t]/2 + vz3[t]/3, -(vz1[t]/3) + vz3[t]/3,
vz1[t]/6 + vz2[t]/2 + vz3[t]/3}
*)


The complete (inhomogeneous) equation with

FD[t_] := f0 Cos[w t]


has the solution

solf = DSolve[eq1, vz[t], t]

(*
Out[38]= {{vz1[t] ->
C[1] Cos[Sqrt[3] t] + C[2] Sin[Sqrt[3] t] - (
f0 (Cos[Sqrt[3] t]^2 Cos[t w] + Cos[t w] Sin[Sqrt[3] t]^2))/(-3 + w^2),
vz2[t] -> C[3] Cos[t] +
C[4] Sin[t] + (-4 xk Cos[t]^2 + 4 w^2 xk Cos[t]^2 +
f0 Cos[t] Cos[t (-1 + w)] + f0 w Cos[t] Cos[t (-1 + w)] +
f0 Cos[t] Cos[t (1 + w)] - f0 w Cos[t] Cos[t (1 + w)] -
4 xk Sin[t]^2 + 4 w^2 xk Sin[t]^2 - f0 Sin[t] Sin[t (-1 + w)] -
f0 w Sin[t] Sin[t (-1 + w)] + f0 Sin[t] Sin[t (1 + w)] -
f0 w Sin[t] Sin[t (1 + w)])/(2 (-1 + w) (1 + w)),
vz3[t] -> C[5] + t C[6] - (f0 Cos[t w])/w^2}}
*)


Here, again, after transforming back to vx, the constants of integration have to be determined by the inial conditions.

I'll leave this last step to the reader.

• Still digesting this; is it correct to assume md = u . ma . Inverse[u]? – Anti Earth Sep 1 '15 at 1:48
• @ Anti Earth md = u.ma.Inverse[u] is not an assumption but a result. Eigensystem[[2]] calculates the Matrix u so that this relation holds, it diagonalizes the Matrix. Exactly that decouples the equations. You might wish to consult reference.wolfram.com/language/tutorial/… – Dr. Wolfgang Hintze Sep 1 '15 at 7:10
• Sorry I rather meant "is it correct to assume that's what the variable md has been assigned?". I've only just noticed its definition now. Why must I solve for the homogeneous and particular solutions separately? Will Mathematica not include the homogeneous solutions automatically in the full solutions? – Anti Earth Sep 1 '15 at 12:32
• @ Anti Earth From your comments I deduce that you should try to get acquainted with the theory of linear differential equations first. Then you'll understand the steps and the aid Mathematica provides as shown. As to you specific question: of course, the full solution includes the homegeneous one. I just have split the two in order to simplify the presentation (maybe not a good idea). – Dr. Wolfgang Hintze Sep 2 '15 at 6:57
• I'm certainly familiar with solving 2nd order constant coefficient in-homogeneous equations, just not a coupled system thereof. I understand why I would solve particular and homogeneous solutions separately, just not why I'd force Mathematica to. – Anti Earth Sep 2 '15 at 8:39