Re: Set of ODE's

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

```
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.