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