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)

 

[alphy]

for your question! It seems like you are on the right track with your approach to solving this system of linear ODEs using Mathematica. Here are a few suggestions to help you solve it using matrix notation:

1. First, make sure that you have defined your matrix X and vector P correctly. This means that the dimensions of X should be NxN (where N is the number of equations) and P should be a vector of length N.

2. Instead of using "solExp == NDSolve[Y'[t] == X.Y[t] && Y[0] = P, X, {t, 0, 10}]", try using "solExp = NDSolve[{Y'[t] == X.Y[t], Y[0] == P}, Y, {t, 0, 10}]". This will return a solution for the vector Y[t] instead of just the matrix X.

3. You can access the individual elements of the vector Y[t] by using Y[t][[1]] for y1(t), Y[t][[2]] for y2(t), and so on. This will make it easier to work with the large set of equations.

4. Lastly, if you are still having trouble with the matrix notation, you can also try using the explicit form of the equations as you mentioned. You can use a loop or a function to define each equation separately, and then use NDSolve to solve the system.

I hope this helps and good luck with your project! Mathematica is a powerful tool for solving systems of ODEs, so with some practice and experimentation, you should be able to successfully solve your problem.