MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Set of ODE's

  • To: mathgroup at
  • Subject: [mg95262] Re: Set of ODE's
  • From: dh <dh at>
  • Date: Tue, 13 Jan 2009 07:03:05 -0500 (EST)
  • References: <gkgovf$3sk$>

Hi Barris,

mathematica evaluates the arguments of a function before calling it. 

That is, it evaluates A.x[t] before calling "+" and "==". From this 

alone, it is not defined if x is a scalar, vector or matrix. Mathematica 

leaves it unevaluated. Next we have A.x + IdentityMatrix[3]. This time 

mathematica  treats A.x[t] obviously as a scalar (I have no clue why).

What can we do? A solution to this is, to declare x as a matrix like e.g.:

x[t_] = Table[ToExpression["x" <> ToString[i] <> ToString[j]][t], {i, 

3}, {j, 3}]

Next problem we face is, that "==" does not thread over matrixes. We 

therefore have to flatten our equations before thtreading. We may define 

an operator for this purpose:

toEq[x_] := x /. {a__List} :> Flatten[{a}] // Thread;

Note further that we need to specify the functions to calculate by 

vars=Flatten@x[t], not only by "x" because this would not be expanded.

Finally we may now write:

A = RandomReal[{1, -1}, {3, 3}];

toEq[x_] := x /. {a__List} :> Flatten[{a}] // Thread;

x[t_] = Table[

   ToExpression["x" <> ToString[i] <> ToString[j]][t], {i, 3}, {j, 3}]

eq = toEq[x'[t] == A.x[t] + IdentityMatrix[3]];

eq = Join[eq, toEq[x[0] == RandomReal[1, {3, 3}]]];

vars = Flatten@x[t];

NDSolve[eq, vars, {t, 0, 20}]

hope this helps, Daniel

Baris Altunkaynak wrote:

> I have a large set of vector/matrix ODE's I want to solve numerically. But I am having some problems with the notation. For example something like that works very well


> A = RandomReal[{1, -1}, {3, 3}];

> NDSolve[{x'[t] == A.x[t], x[0] == RandomReal[1, {3, 3}]}, x, {t, 0, 20}]


> but it doesn't work if I add a constant term to the right hand side of the equation:


> NDSolve[{x'[t] == A.x[t] + IdentityMatrix[3], x[0] == RandomReal[1, {3, 3}]}, x, {t, 0, 20}]


> I guess the problem is that Mathematica adds A.x[t] to all the components of the IdentityMatrix[3] as if it is a scalar quantity. I defined the identity matrix to be a matrix function with zero derivative etc, so that it doesn't get evaluated immediately and that works but there should be a better way of doing it with the Hold function I guess. Do you have any idea how to do this in a more elegant way?


> thank you,

> Baris Altunkaynak





  • Prev by Date: Re: floating point comparison
  • Next by Date: NEWBIE QUESTION - How to create a barchart with two data sets.
  • Previous by thread: Re: How to stop the latest Mathematica hijacking the .nb extension in
  • Next by thread: NEWBIE QUESTION - How to create a barchart with two data sets.