NDSolve in matrix system with Mathematica

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

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

with a "pseudo" Dirac Delta initial condition:

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

This system can be solved with NDSolve:

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:

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

xi[1, 0.1] /. sol

but the answer was simply

{xi[1, 0.1]}

. Using Flatten I obtain the same result (without the braces) and if I try with

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.


Arturo.

 

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

n = 4;

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

X0 = Table[0, {n}];
X0[[Round[n/2]]] = 1;

e.g.

M = {{1, 2, 0, 0}, {0, 2, 3, 0}, {0, 0, 1, 2}, {1, 0, 0, 1}};

diffEq = Table[X'[t][[i]] == -(M.X[t])[[i]], {i, 1, n}];
initialCond = Table[X[0][[i]] == X0[[i]], {i, 1, n}];

sol = NDSolve[{Join[diffEq, initialCond]}, X[t], {t, 0, 1}]

And you can plot solutions with:

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

 

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

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

xi[2,0.1]

or

{xi[2,0.1]}

I already took a look to

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.

 

[djelovin]

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

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

In my example

x[2, 0.1]

returns

0.818774

 

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