Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

NDSolve in matrix system with Mathematica

  1. Mar 28, 2012 #1
    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):
    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,
  2. jcsd
  3. Mar 28, 2012 #2
    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:

    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;
    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
  4. Mar 29, 2012 #3
    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):

    but I always obtain
    Code (Text):

    Code (Text):

    I already took a look to
    Code (Text):

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

  5. Mar 29, 2012 #4
    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]
  6. Mar 30, 2012 #5
    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.

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook