matrix diffential equations follow up
- To: mathgroup at smc.vnet.net
- Subject: [mg27727] matrix diffential equations follow up
- From: "Martin Richter" <mrMICE.fi at cbs.dk>
- Date: Tue, 13 Mar 2001 03:52:45 -0500 (EST)
- Organization: UNI-C
- References: <MDzp6.2242$96.11689@ralph.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi
Thanks to anyone who did come up with suggestions for defining matrix
differential equations. I haven't have the time to look at it before now.
I think it was Marc Harder who come up with the solution, I just have to
change one line. I hav just paste my code, it straight forward to put in,
if anyone want to do it.
What I wanted to achieve was to solve the matrix equation
dPhi/dt = L(t)dt, L(s)=I (identity)
If L(t)==L the solution is Exp((t-s)*L)
If L(t) and L(s) commutes one get Exp(\int_0^t L(u)du)
But if the matrices doesn't commute (at all time) there is not in general a
close form solution. My case is such a case,
so my task was simply to compare the error from different approximations
methods of the equation. It turn out to be
a more complicated programming than expected (for a newbie in Mathematica),
and the problem didn't have a high priority for the project.
Martin
<< LinearAlgebra`MatrixManipulation`;
Clear[\[Alpha], \[Kappa], \[Lambda]P, \[Lambda]\[Delta], \[Lambda]v,
\[Sigma]\
\[Delta], \[Sigma]v, \[Nu], \[Nu]11, \[Nu]12, \[Nu]21, \[Nu]22]
\[Alpha] = 1.8952; \[Kappa] = 3.4648;
\[Lambda]P = 0.1108; \[Lambda]\[Delta] = -0.0233; \[Lambda]v = -0.0651;
\[Sigma]\[Delta] = 1.8387; \[Sigma]v = 0.6126;
\[Nu]11 = -0.5915; \[Nu]12 = -0.2511; \[Nu]21 = 0.1737; \[Nu]22 = 0.1880;
Clear[L, \[CapitalPhi], u]
\[Nu][t_] := \[Nu]11*Cos[2*Pi*t] + \[Nu]12*Sin[2*Pi*t] + \[Nu]21*
Cos[4*Pi*t] + \[Nu]22*Sin[4*Pi*t];
L[t_] := {{0, -1, -1/2 + \[Lambda]P*
Exp[1/2*\[Nu][
t]]}, {0, -\[Alpha], \[Lambda]\[Delta]*\[Sigma]\[Delta]*
Exp[1/2*\[Nu][t]]}, {0, 0, -\[Kappa] + \[Lambda]v*\[Sigma]v}};
fSymbol[f_String, i_Integer, j_Integer] :=
ToExpression[f <> ToString[i] <> ToString[j]];
\[CapitalPhi][t_] :=
Table[fSymbol["\[CapitalPhi]", i, j][t], {i, 3}, {j, 3}];
(* Now come the extra line*)
\[CapitalPhi]Symbol := Table[fSymbol["\[CapitalPhi]", i, j], {i, 3}, {j,
3}];
eqns = Flatten[{MapThread[Equal, {\[CapitalPhi]'[u], L[u].\[CapitalPhi][u]},
2], MapThread[Equal, {\[CapitalPhi][s], IdentityMatrix[3]}, 2]}];
Clear[s, t]
s = 1; (* starting point*)
t = 1.1; (* end point *)
(* a small change in the next line *)
Sol = NDSolve[eqns, Flatten[\[CapitalPhi]Symbol], {u, s, t}];
Plot[Evaluate[\[CapitalPhi]22[u] /. Sol], {u, s, t}];
Plot[Evaluate[MatrixNorm[Flatten[\[CapitalPhi][u]] /. Sol, 2]], {u, s, t},
PlotRange -> {1, 2}];