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

Question

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!

Answer 1

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.