NDSolve and matrix formulation do not produce the same result. Why and which is more precise?
It seems you have made some error, the eigensystem approach is shown to work:
First, You may find it useful to know that NDSolve
can work directly with vector equations: (This is also validating the first part of your manipulation)
a={{8, 0, -8}, {0, 12, -4}, {-8, -4, 20}} 10^5
b={{1, 1, 0}, {1, 1, 0}, {0, 0, 0}} 1770.17
c={{-32, -22, 10}, {-22, -20, 10}, {10, 10, -8}} 7.57418
f={-3, -2, 1} 1415.29
ia = Inverse[a];
ic = Inverse[c];
sol = First@NDSolve[ {
u''[t] == (ic.a).u[t] + (ic.b).u'[t] ,
u[0] == ia.f,
u'[0] == {0, 0, 0}
}, u[t], {t, 0, 5}]
aa[t_, i_] := (u[t] /. sol)[[i]] - (ia.f)[[i]]
Plot[aa[t, 1], {t, 0, 1.5}, PlotRange -> All]
your first order form gives the same result:
h = {{ConstantArray[0, {3, 3}], IdentityMatrix[3]}, {ic.a, ic.b}};
sol = First@NDSolve[{
{u'[t], v'[t]} == {
h[[1, 1]].u[t] + h[[1, 2]].v[t],
h[[2, 1]].u[t] + h[[2, 2]].v[t]},
{u[0], v[0]} == {ia.f, {0, 0, 0}}}, {u[t], v[t]}, {t, 0, 5}]
I think this capability is fairly new.. (this is v10.1). unfortunately it doesn't appear to generalize to higher order tensor this throws an error:
NDSolve[{z'[t] == h.z[t], z[0] == {ia.f, {0, 0, 0}}}, z[t], {t, 0, 5}]
we can however pose this as a 6-vector, using a 6x6 form of h:
h66 = {
Join[h[[1, 1, 1]], h[[1, 2, 1]]],
Join[h[[1, 1, 2]], h[[1, 2, 2]]],
Join[h[[1, 1, 3]], h[[1, 2, 3]]],
Join[h[[2, 1, 1]], h[[2, 2, 1]]],
Join[h[[2, 1, 2]], h[[2, 2, 2]]],
Join[h[[2, 1, 3]], h[[2, 2, 3]]]};
(* surely some incantation of Flatten[h,levelspec] does the same..*)
sol = First@
NDSolve[{z'[t] == h66.z[t], z[0] == Flatten[{ia.f, {0, 0, 0}}]},
z[t], {t, 0, 5}]
aa[t_, i_] := (z[t] /. sol)[[i]] - (ia.f)[[i]]
Plot[aa[t, 1], {t, 0, 1.5}, PlotRange -> All]
same plot
Now we can apply your eigensystem solution:
{eval, evec} = Eigensystem[h66];
csol = First@
Solve[ Sum[ cc[k] evec[[k]] , {k, 6}] ==
Flatten[{ia.f, {0, 0, 0}}] , Array[cc, 6]];
eigsol = Sum[ (cc[k] /. csol) evec[[k]] Exp[eval[[k]] t], {k, 6}];
Plot[Table[Re[eigsol[[k]]] - (ia.f)[[k]], {k, 3}], {t, 0, 2},
PlotRange -> All]
Exact same plot ( no noise ! )
limit value Chop[eigsol[[1]] - (ia.f)[[1]] /. t -> Infinity]
0.00884556
which is First@(-Inverse[a].f)
So Mathematica can, in fact, solve your equations exactly using DSolve
instead of NDSolve
:
solution = DSolve[equations, variables, t];
koti = {\[CurlyPhi]1[t], \[CurlyPhi]2[t], \[CurlyPhi]3[t]} /. solution;
Plot[Evaluate[Re /@ koti], {t, 0, 1}, PlotRange -> Full,
AxesLabel -> {t, None}, PlotPoints -> 100, ImageSize -> Large]
Note that I removed the Rationalize
from equations
before feeding it into DSolve
. If you don't, Mathematica will try to work with exact roots of sixth-order polynomials, and that way lies madness.
What I noticed in this graph is that the green graph is not as smooth as one would expect. This non-smoothness persists even if you increase the number of plot points (note the high value of PlotPoints
in the above code) and if you zoom in on the graph:
This looks like a large oscillation with a very small oscillation superposed on it. Notably, in your second method, eigenVal
works out to be
(* {-11.2869 + 344.474 I, -11.2869 - 344.474 I, -9.99201*10^-15 + 162.498 I,
-9.99201*10^-15 - 162.498 I, -5.40686 + 40.5931 I, -5.40686 - 40.5931 I} *)
and $2\pi/344.474 \approx 0.018$, which matches up nicely with the period of the small oscillations you can see on the graph above.
What I suspect is happening here is that your system has a high-frequency mode whose corresponding amplitude $C_k$ is low in the "true" solution. However, something is going wrong with your application of boundary conditions; resulting in this high-frequency oscillation being assigned a much higher amplitude $C_k$, making it noticeable.
In conclusion: the results of NDSolve
and DSolve
are generally consistent with each other. Moreover, a close examination of those two results has convinced me that there is an error in your home-brew method (in particular, the application of the initial conditions) that results in the weird-looking solution that is generated.
As an aside: DSolve
is probably using a version of your algorithm to solve the equation exactly. But when I have a choice between using pre-tested algorithms to accomplish a task versus writing my own, I'll use the pre-tested algorithms every time.