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