Matrix Systems of ODEs - Mathematica

[brydustin]

Matrix Systems of ODEs -- Mathematica

I'm trying to solve the classic "systemm of linear ODEs" of the form: Y(t)' = X*Y(t)
Its homogenous, so it wouldn't hurt to rewrite it as Y(t)' - X*Y(t) = 0 (if that helps?)


so here is my attempt at the solution

solExp == NDSolve[Y'[t] == X.Y[t] && Y[0] = P, X, {t, 0, 10}]

where P is a vector of numbers (they are each fixed) -- its my list of initial conditions, because each equation is first order, each eq-n needs one initial condition.

you can imagine the situation of:

(dy1/dt)= x1,1*y1(t) + x1,2*y2(t) + ... x1,N*yN(t)
(dy2/dt)= x2,1*y1(t) + x2,2*y2(t) + ... x2,N*yN(t)
.
.
.
(dyN/dt) = xN,1*y1(t) + xN,2*y2(t) + ... xN,N*yN(t)

we are given the intial condition that the set Y[t=0] = P = {y1(0), y2(0),...yN(0)}
I can't seem to get this to work using matrix notation... although I know how to do it explicitly for a few equations, where you list each on in the code... the thing is: I am dealing with a LARGE set of equations and so this wouldn't be practical.

Thanks

 

[Bill Simpson]

Suggestion: First figure out all the errors you have in creating each of the arguments you are going to pass to NDSolve before you even consider using NDSolve.

In other words practice on
eqns={something}
until you get something exactly right and the result exactly matches the syntax of
http://reference.wolfram.com/mathematica/ref/NDSolve.html
for systems of equations.

Then practice on
vars={somethingelse}
until you get somethingelse exactly right and the result exactly matches the syntax of
http://reference.wolfram.com/mathematica/ref/NDSolve.html
for systems of equations.

Then you might be ready for NDSolve[eqns,vars,{t,0,10}]

 

[Simon_Tyler]

Hi brydustin,

Bill's right... there's a lot wrong in your original post.

Here's a simple example of a 2*2 constant coefficient system to help get you started:

X = Array[x, {2, 2}];      (* Arbitrary constant coefficient matrix *)
Y[t_] := {y1[t], y2[t]};   (* functions to solve for *)
Y0 = {y10, y20};           (* initial conditions  *)

(*  Generate equations *)
eqns = Thread /@ {Y'[t] == X.Y[t], Y[0] == Y0}
(* Alternately use something like *)
(* eqns = And@@Flatten[Thread /@ {Y'[t] == X.Y[t], Y[0] == Y0}]  *)
(* Or *)
(* eqns = Thread[Flatten[{Y'[t] - X.Y[t], Y[0] - Y0}] == 0] *)(* Solve -- will produce a mess  *)
soln = DSolve[eqns, Y[t], t];

(*  make the output a little nicer on th eyes...  *)
soln /. x[1, 1]^2 + 4 x[1, 2] x[2, 1] - 2 x[1, 1] x[2, 2] + 
    x[2, 2]^2 -> \[Delta]^2 // FullSimplify[#, \[Delta] > 0] &

Also, if it really is a system of constant coefficient 1st order DEs, then it might be better to use MatrixExp to generate the solution, this way you can minimize the numerical errors that are inevitable in the NDSolve: The following will return "True"

MatrixExp[t X].Y0 == (Y[t] /. First[soln]) // Simplify

 

[brydustin]

Hi brydustin,

Bill's right... there's a lot wrong in your original post.

Here's a simple example of a 2*2 constant coefficient system to help get you started:

X = Array[x, {2, 2}];      (* Arbitrary constant coefficient matrix *)
Y[t_] := {y1[t], y2[t]};   (* functions to solve for *)
Y0 = {y10, y20};           (* initial conditions  *)

(*  Generate equations *)
eqns = Thread /@ {Y'[t] == X.Y[t], Y[0] == Y0}
(* Alternately use something like *)
(* eqns = And@@Flatten[Thread /@ {Y'[t] == X.Y[t], Y[0] == Y0}]  *)
(* Or *)
(* eqns = Thread[Flatten[{Y'[t] - X.Y[t], Y[0] - Y0}] == 0] *)


(* Solve -- will produce a mess  *)
soln = DSolve[eqns, Y[t], t];

(*  make the output a little nicer on th eyes...  *)
soln /. x[1, 1]^2 + 4 x[1, 2] x[2, 1] - 2 x[1, 1] x[2, 2] + 
    x[2, 2]^2 -> \[Delta]^2 // FullSimplify[#, \[Delta] > 0] &

Also, if it really is a system of constant coefficient 1st order DEs, then it might be better to use MatrixExp to generate the solution, this way you can minimize the numerical errors that are inevitable in the NDSolve: The following will return "True"

MatrixExp[t X].Y0 == (Y[t] /. First[soln]) // Simplify

Cool! This worked, but I also tried the following:
d = Dimensions[W][[1]];
GeneralSolution = 0;
Eigen = Eigensystem[W];
GeneralSolution = 0;
m = Transpose[Last[Eigen]];
CValues = LinearSolve[Transpose[Last[Eigen]], P]; (*P is the intial condition vector*)
For[i = 0, i < d,
GeneralSolution =
GeneralSolution + CValues[]*E^(Eigen[[1, i]]*t)*Eigen[[2, i]], i++]

Which seems to work, and is equivalent to Sum(c*v*e^(\(lambda)*t)), where c's are the constant coefficients for a particular sol-n, {v,\lamdba} is the corresponding eigenpair, and t is time. I'm pretty sure this was the way we learned in intro ode, which is equivalent to the matrix exponential (which I did before, but crashed sometime, for some reason)