Let's assume a vector of variables $\vec a=(a_1(t),a_2(t),a_3(t))$ and let's assume a following DE$$C\ddot{\vec a}=B\dot{\vec a}+A\vec{a}+\vec F.$$ If $A$ and $C$ are both reversible, than the DE can be written as $$\ddot{\vec a}=C^{-1}B\dot{\vec a}+C^{-1}A(\vec a+A^{-1}\vec F)$$ now making a crucial step $${(\vec a+A^{-1}\vec F)}''=C^{-1}B{(\vec a+A^{-1}\vec F)}'+C^{-1}A(\vec a+A^{-1}\vec F)$$ where $'$ still denotes a time derivative - it's just a cleaner syntax. By saying $(\vec a+A^{-1}\vec F)=\vec u$ we derived a very compact form of the second ordered DE system $${\vec u }''=C^{-1}B{\vec u}'+C^{-1}A\vec u.$$ To reduce the order from second to first we define yet another vector $\vec v={\vec u}'$. Doing so, our system can be rewritten in the following form $${\binom{\vec u}{\vec v}}'=\begin{pmatrix} 0 &I \\ C^{-1}A&C^{-1}B \end{pmatrix}\binom{\vec u}{\vec v}$$ which is in a $${\vec z}'=H\vec z$$ form and can be solved via eigensystems theories $$\vec z (t)=\sum C_k\vec \lambda _ke^{\lambda_kt}$$ for eigenvalues $\lambda_k$ and eigenvectors $\vec \lambda_k$.
PROBLEM
The system can be "easily" solved via NDSolve[]. The equations are
equations = {800000. φ1[t] - 800000. φ3[t] +
0.105 (16858.8 Derivative[1][φ1][t] +
16858.8 Derivative[1][φ2][t] +
229 (-176.58 + 10.08 φ1''[t] +
6.93 φ2''[t] - 3.15 φ3'' [t])) == 0,
1.2*10^6 φ2[t] - 400000. φ3[t] +
0.105 (16858.8 Derivative[1][φ1][t] +
16858.8 Derivative[1][φ2][t] +
229 (-117.72 + 6.93 φ1''[t] +
6.3 φ2''[t] - 3.15 φ3''[t])) == 0,
1415.29 - 800000. φ1[t] - 400000. φ2[t] +
2.*10^6 φ3[t] - 75.7418 φ1''[t] -
75.7418 φ2''[t] + 60.5934 φ3''[t] == 0};
and using "brutal" NDSolve[]
n = 4;
variables =
Table[ToExpression["φ" <> ToString[i]], {i, 1, n}];
For[i = 1, i <= n - 1, i++,
AppendTo[equations, variables[[i]][t] == 0 /. t -> 0]];
For[i = 1, i <= n - 1, i++,
AppendTo[equations, D[variables[[i]][t], t] == 0 /. t -> 0]];
solution =
Quiet[NDSolve[Rationalize[equations],
Table[variables[[i]][t], {i, 1, n - 1, 1}], {t, 0, 5}]];
koti = Table[solution[[1, i, 2]], {i, 1, n - 1, 1}];
Plot[koti, {t, 0, 1.5}, PlotRange -> All, PlotLegends -> Automatic,
AxesLabel -> {"t [s]", "φ [rad]"},
BaseStyle -> {FontFamily -> "Courier New", FontSize -> 10}]
produces
HOWEVER using all the theory written above and the code below produces something... well something that does not look the same. Make sure you clear the "equations" because at this point boundary conditions were added.
matrixSecondD =
Normal@CoefficientArrays[
equations, {φ1''[t], φ2''[t], φ3''[
t]}][[2]];
matrixFirstD = -Normal@
CoefficientArrays[
Normal@CoefficientArrays[
equations, {φ1''[t], φ2''[
t], φ3''[t]}][[1]], {φ1'[
t], φ2'[t], φ3'[t]}][[2]];
matrixZeroD = -Normal@
CoefficientArrays[
Normal@CoefficientArrays[
Normal@
CoefficientArrays[
equations, {φ1''[t], φ2''[
t], φ3''[t]}][[1]], {φ1'[
t], φ2'[t], φ3'[t]}][[
1]], {φ1[t], φ2[t], φ3[t]}][[2]];
matrixNonhomogeneousPart = -Normal@
CoefficientArrays[
Normal@CoefficientArrays[
Normal@CoefficientArrays[
equations, {φ1''[t], φ2''[
t], φ3''[t]}][[1]], {φ1'[
t], φ2'[t], φ3'[t]}][[
1]], {φ1[t], φ2[t], φ3[t]}][[1]];
matrixSecondD.{φ1''[t], φ2''[t], φ3''[
t]} == matrixFirstD.{φ1'[t], φ2'[
t], φ3'[t]} +
matrixZeroD.{φ1[t], φ2[t], φ3[t]} +
matrixNonhomogeneousPart;
id = IdentityMatrix[3];
zeroes = ConstantArray[0, {3, 3}];
matrix = ArrayFlatten[{{zeroes,
id}, {Inverse[matrixSecondD].matrixZeroD,
Inverse[matrixSecondD].matrixFirstD}}];
eigenVal = Eigensystem[matrix][[1]];
eigenVec = Eigensystem[matrix][[2]];
constants = Table[ToExpression["C" <> ToString[i]], {i, 1, 2 (n - 1)}];
freeMatrix = Inverse[matrixZeroD].matrixNonhomogeneousPart;
equationSolution =
Sum[Table[
constants[[i]]*Exp[eigenVal[[i]]*t] eigenVec[[i]] -
If[i <= (n - 1), freeMatrix[[i]], 0], {i, 1, Length[eigenVal],
1}][[j]], {j, 1, 2 (n - 1)}];
equationsForConst =
Table[equationSolution[[i]] == 0 /. t -> 0, {i, 1,
Length[eigenVal]}];
system = Solve[equationsForConst, constants];
constants = constants /. Flatten[system];
equationSolution =
equationSolution /.
Thread[Table[
ToExpression["C" <> ToString[i]], {i, 1, 2 (n - 1)}] ->
constants];
koti = Drop[equationSolution, -(n - 1)];
Plot[Evaluate[ReIm /@ koti], {t, 0, 1}, PlotRange -> Full,
AxesLabel -> {t, None},
PlotStyle -> (Sequence @@ {Directive[Thin, ColorData[1][#]],
Directive[Dashed, ColorData[1][#]]} & /@ Range[Length@koti]),
PlotLegends ->
Column@{LineLegend@*Sequence @@
Transpose[{ColorData[1][#], "Eq. " <> ToString@#} & /@
Range[Length@koti]],
LineLegend[{Thin, Dashed}, {"Real", "Imaginary"}]},
ImageSize -> Large]
Why are solutions to the DE system not the same in both cases? And which one is the correct one? NOTE that this is the first time I am working with matrices in mathematica, so there is probably a really stupid mistake somewhere.
ps: that last part (plotting the real and imaginary party) is a copy paste code from Edmund in THIS answer.
a = (a+A^(-1) F)Am I missing something? $\endgroup$(φ1^′′)[t], which is definitely not right. $\endgroup$u = (a+A^(-1) F)and not something else? Because with this definition I don't change the equation. Note thatA^(-1) Fis a constant vector. @MichaelSeifert: That's a notation this forum uses. If you copy paste the code in Mathematica, it should give you the derivatives. $\endgroup$