# NDSolve in matrix system with Mathematica

1. Mar 28, 2012

### jatorre

Hello everyone.

This is my first post after looking for information about a simple question. Unfortunately I didn't found information in any other posts in this forum so, here we go!

I'm trying to solve numerically a set of N ODEs with Mathematica 7.0. Let me call x(t) an array of N elements and M a NxN matrix. I would like to solve the problem dx/dt = - M x.

So I firstly define the array:

Code (Text):

X[t_] = Table[xi[i, t], {i, N}];

with a "pseudo" Dirac Delta initial condition:

Code (Text):

X0 = SparseArray[Table[{i} -> 0, {i, N}]];
X0[[N/2]] = 1;

This system can be solved with NDSolve:

Code (Text):

eqns = Thread /@ {X'[t] == -M.X[t], X[0] == X0};
sol = NDSolve[eqns, X[t], {t, 0, 1}];

If I plot the solution I can check the correct answer:

Code (Text):

Plot[Evaluate[Table[xi[i, t] /. sol, {i, N}]], {t, 0, 1}]

And here comes my question. How to obtain a specific value of xi (e.g., xi[1,0.1])?

I tried to do it using directly
Code (Text):
xi[1, 0.1] /. sol
but the answer was simply
Code (Text):
{xi[1, 0.1]}
. Using Flatten I obtain the same result (without the braces) and if I try with
Code (Text):
X[0.1]/.sol
I obtain an array of elements xi[i,0.1] with i=1,N.

At this point I have no idea of how to deal with this. I am quite sure that this has to be a silly question but, as far as I know, I cannot obtain the correct answer.

Thank you very much in advanced and my apologies if this is a duplicated question.

Best regards,
Arturo.

2. Mar 28, 2012

### djelovin

Your syntax is wrong in many ways and you can't use N as a variable name since it's Mathematica's standard function. Here is your code correctly rewrited:

e.g.
Code (Text):
n = 4;
Code (Text):
X[t_] = Table[xi[i, t], {i, 1, n}];
Code (Text):
X0 = Table[0, {n}];
X0[[Round[n/2]]] = 1;
e.g.
Code (Text):
M = {{1, 2, 0, 0}, {0, 2, 3, 0}, {0, 0, 1, 2}, {1, 0, 0, 1}};
Code (Text):
diffEq = Table[X'[t][[i]] == -(M.X[t])[[i]], {i, 1, n}];
initialCond = Table[X[0][[i]] == X0[[i]], {i, 1, n}];
Code (Text):

sol = NDSolve[{Join[diffEq, initialCond]}, X[t], {t, 0, 1}]
And you can plot solutions with:
Code (Text):
Plot[Evaluate[Table[xi[i, t] /. sol, {i, 1, n}]], {t, 0, 1}]

Last edited: Mar 28, 2012
3. Mar 29, 2012

### jatorre

Dear djelovin,

thank you very much for the quick answer. I chose very bad names for the example. The syntax you propose is indeed much clear than mine, thanks :)

Anyway, I still cannot print a value for xi (eg., xi[2,0.1]). The plot is already done and correctly shown, but if I try to print a value Mathematica always returns the input. I tried:

Code (Text):

xi[2,0.1]/.sol
Evaluate[xi[2,0.1]/.sol
Evaluate[xi[2,0.1]/.Flatten[sol]
Evaluate[X[0.1][[2]]/.sol
Evaluate[X[0.1][[2]]/.sol[[2]]

but I always obtain
Code (Text):

xi[2,0.1]

or
Code (Text):

{xi[2,0.1]}

I already took a look to
Code (Text):

InputForm[sol]

but all I obtained was what seems to be a nested array in which I am not able to recover any specific data.

Thanks again for your attention :)

Regards,
Arturo.

4. Mar 29, 2012

### djelovin

You have to use mathematica rules for t variable too. If you wont function that returns numerical value, you can define if like this:

Code (Text):
x[2, x_] := (xi[2, t] /. sol /. t -> x)[[1]]

In my example
Code (Text):
x[2, 0.1]
returns

5. Mar 30, 2012

### jatorre

Dear djelovin,

Thank you for the clarification. Now it works perfectly :) After reading the documentation about rules, I now understand how they work and how to obtain (as in your example) numerical values.

Regards,
Arturo.