Results 1 to 9 of 9

Math Help - Mathematica - ODE system solving problem

  1. #1
    Newbie
    Joined
    Nov 2009
    Posts
    8

    Question Mathematica - ODE system solving problem

    Hello All,

    Using Mathematica I met a problem with solving of ODE. Maybe this problem has a non-technical, mathematical character

    First step I needed to solve this system (kinetic equations for isotopic exchange)

    /parameters calculated from the experimental conditions/
    K1=0.107
    K2=0.541
    K3=0.844
    α0=0.981
    αs0=0.0024

    result := DSolve[{α'[t] == (α0 - α[t])*K1 + (αs[t] - α[t])*K2, αs'[t] == (α[t] - αs[t])*K3, α[0] == αs0, αs[0] == αs0, {α[t], αs[t]}, t]

    α00[t_] = α[t]/.result

    It has an analytical solution - sum of exponents, and α0[t] behaves itself as functions like f(t)=1-p*exp(-qt), 0<p<1, q>0 on interval {0,+inf}. Also I found numerical solution using NDSolve procedure, interpolating function on the required interval {0,100} looked like the same as analytical one.

    Then I needed to solve the similar system:

    result := NDSolve[{α'[t] == (α00[t] - α[t])*K1 + (αs[t] - α[t])*K2, αs'[t] == (α[t] - αs[t])*K3, α[0] == αs0, αs[0] == αs0, {α[t], αs[t]}, {t,0,100}]

    The difference from the first system is contained in replacing of the constant "α0" by the function "α00[t]", which is a solution of the first system.

    The solution of the second system for α[t] must look like the aforesaid exponent (I know it because I saw the solution in other software environment), but in the output there is written "The derivative is not consistent with initial conditions" or something like that, and i cannot plot α[t] (obtained from the second system solution) curve. I don't know what is wrong, I'm not an expert in math. I tried to put "α00[t]" as both analytical and numerical solution, but result was the same.

    I saw the "Help", but didn't find any solution for my case.

    Could you give me your suggestions how to solve this system? Maybe I need to change some parameters?

    I work in Mathematica 5.2 for students

    Thank you for your help in advance!
    Last edited by LeChat; November 25th 2009 at 02:43 PM.
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Hi. I think the \alpha s function in the first DSolve is conflicting with the \alpha s in the NDSolve. So I changed the first to \alpha 2s. I changed a few other things too. Try and work through my code and see what I did. Not sure why the alpha is not converting over when I convert it to raw form to cut and paste. You'll probably have to manually convert it or just change the function names manually in your code. Also, don't start user-defined variable names with upper case. They may conflict with internal commands. Personally, the alpha thing is too hard to input for me, I would just as well have called the functions u(t), v(t), w(t), s(t) and so forth to make it easier to track in Mathematica and avoid naming conflicts. Your call. Anyway, I'm using ver. 7 and the code below worked fine. Hope it does so with 5.4.

    Code:
    K1 = 0.107
    K2 = 0.541
    K3 = 0.844
    \[Alpha]0 = 0.981
    \[Alpha]s20 = 0.0024
    mySolutions = DSolve[
       {Derivative[1][\[Alpha]][t] == 
         (\[Alpha]0 - \[Alpha][t])*K1 + 
          (\[Alpha]s2[t] - \[Alpha][t])*K2, 
        Derivative[1][\[Alpha]s2][t] == 
         (\[Alpha][t] - \[Alpha]s2[t])*K3, 
        \[Alpha][0] == \[Alpha]s20, \[Alpha]s2[0] == 
         \[Alpha]s20}, {\[Alpha], \[Alpha]s2}, t]
    Plot[{\[Alpha][t], \[Alpha]s2[t]} /. mySolutions, 
      {t, 0, 10}]
    myf[t_] = First[Evaluate[\[Alpha]s2[t] /. 
         mySolutions]]
    sol = NDSolve[{Derivative[1][\[Alpha]][t] == 
         myf[t]*K1 + (\[Alpha]s[t] - \[Alpha][t])*K2, 
        Derivative[1][\[Alpha]s][t] == 
         (\[Alpha][t] - \[Alpha]s[t])*K3, 
        \[Alpha][0] == \[Alpha]s20, \[Alpha]s[0] == \[Alpha]s20}, 
       {\[Alpha], \[Alpha]s}, {t, 0, 10}]
    Plot[{\[Alpha][t], \[Alpha]s[t]} /. sol, {t, 0, 10}]
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Nov 2009
    Posts
    8
    Ok! Good help, shawsend, it works! Thanks a lot! And thanks especially for the fast answer.
    I took into account all the remarks:
    Code:
    k1=0.844
    k2=0.541
    k3=0.107
    k4=0.981
    k5=0.0024
    k6=0.0024
    
    result := NDSolve[{u'[t] == (k4 - u[t])*k1 + (w[t] - u[t])*k2, w'[t] == (u[t] - w[t])*k3, u[0] == k5, w[0] == k5}, {u[t], w[t]}, {t, 0, k6}]
    
    u01[t_] = First[Evaluate[u[t] /. result]]
    
    result := NDSolve[{u'[t] == (u01[t] - u[t])*k1 + (w[t] - u[t])*k2, w'[t] == (u[t] - w[t])*k3, u[0] == k5, w[0] == k5}, {u[t], w[t]}, {t, 0, k6}]
    
    u02[t_] = First[Evaluate[u[t] /. result]]
    Next question: I need to repeat this procedure several times, but i have certain dificulties with function "For":

    Code:
    For[i = 2; uw[t] = u01[t], i < 4, i++, result := NDSolve[{u'[t] == (uw[t] - u[t])*k1 + (w[t] - u[t])*k2, w'[t] == (u[t] - w[t])*k3, u[0] == k5, w[0] == k5}, {u[t], w[t]}, {t, 0, k6}]; uw[t] = First[Evaluate[u[t] /. result]]]
    This code doesn't evaluate function uw[t]. (checked by uw[10])
    In what place I was wrong here?
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Super Member
    Joined
    Aug 2008
    Posts
    903
    The best way to enumerate in Mathematica is with the Table construct. In the code below, I encapsulate the two NDSolves in a Table construct called mySolutions. Note I need semicolons after each command for this. I just varied k1 in the first NDSolve with the command k1*n just for fun. The table returns a 2D array, each row has the two solutions for each value of n which I can enumerate with mySolutions[[n,m]] or in the case below, I plotted all the solutions by passing mySolutions to Evaluate in the final plot command. Notice the final {u01[t], u02[t]} which is the entry actually stored in the table.

    You can cut and paste Mathematica code as text as long as you first convert the cell to raw format via select the cell, then choose Cell/Convert To/Raw. It strips all the special format but then you can move it in and out of text displays.

    Code:
    mySolutions = Table[
       result := NDSolve[
          {Derivative[1][u][t] == 
            (k4 - u[t])*(k1*n) + 
             (w[t] - u[t])*k2, 
           Derivative[1][w][t] == 
            (u[t] - w[t])*k3, u[0] == k5, 
           w[0] == k5}, {u[t], w[t]}, 
          {t, 0, k6}]; u01[t_] = 
         First[Evaluate[u[t] /. result]]; 
        result := NDSolve[
          {Derivative[1][u][t] == 
            (u01[t] - u[t])*k1 + 
             (w[t] - u[t])*k2, 
           Derivative[1][w][t] == 
            (u[t] - w[t])*k3, u[0] == k5, 
           w[0] == k5}, {u[t], w[t]}, 
          {t, 0, k6}]; u02[t_] = 
         First[Evaluate[u[t] /. result]]; 
        {u01[t], u02[t]}, {n, 1, 5}]
    Plot[Evaluate[mySolutions], {t, 0, k6}]
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Newbie
    Joined
    Nov 2009
    Posts
    8
    I'm sorry, it seems I didn't understand Your last offer (of course, it works!) or maybe You didn't understand what do I want. I wanted to use a recursion in a solution, when solution of the first system is placed to the next system, and to repeat it for n times:

    Code:
    result := NDSolve[
          {Derivative[1][u][t] == 
            (k4 - u[t])*k1 + 
             (w[t] - u[t])*k2, 
           Derivative[1][w][t] == 
            (u[t] - w[t])*k3, u[0] == k5, 
           w[0] == k5}, {u[t], w[t]}, 
          {t, 0, k6}]; u01[t_] = 
         First[Evaluate[u[t] /. result]]; 
        result := NDSolve[
          {Derivative[1][u][t] == 
            (u01[t] - u[t])*k1 + 
             (w[t] - u[t])*k2, 
           Derivative[1][w][t] == 
            (u[t] - w[t])*k3, u[0] == k5, 
           w[0] == k5}, {u[t], w[t]}, 
          {t, 0, k6}]; u02[t_] = 
         First[Evaluate[u[t] /. result]]
        result := NDSolve[
          {Derivative[1][u][t] == 
            (u02[t] - u[t])*k1 + 
             (w[t] - u[t])*k2,
           Derivative[1][w][t] == 
            (u[t] - w[t])*k3, u[0] == k5, 
           w[0] == k5}, {u[t], w[t]}, 
          {t, 0, k6}]; u03[t_] = 
         First[Evaluate[u[t] /. result]]
        ...       
        result := NDSolve[
          {Derivative[1][u][t] == 
            (u0/n-1/[t] - u[t])*k1 + 
             (w[t] - u[t])*k2,
           Derivative[1][w][t] == 
            (u[t] - w[t])*k3, u[0] == k5, 
           w[0] == k5}, {u[t], w[t]}, 
          {t, 0, k6}]; u0/n/[t_] = 
         First[Evaluate[u[t] /. result]]
    And at the end only u0/n/[t] will be plotted.

    I don't see how can i do it using Table construct. How can i write a link on the previous cell in one-dimensional table inside this table, if it's possible?

    Or maybe it would be better to organize this procedure in cycle such as "For", or "While"? But all my attempts to do it failed
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Super Member
    Joined
    Aug 2008
    Posts
    903
    Ok, so you're passing the results from one NDSolve, into the next NDSolve. This is what may work but really, if it were mine, I'd first test the procedure with a simple setup that I know what the answer, make sure it gives me that answer, then run this setup. So I get the first solutions and call it myCurrentFunction. Then I set up a table and run NDSolve 8 times, and the end of each one, I set myCurrentFunction to the current solution of NDSolve so that it can be passed to the next iteration. I save all the solutions in myFunctionTable than prepend the first one to the table, then iterate the plots of all 9 with the final plot command and display the set of 9 in a graphics array. Not sure why those 6 and more are flat lining though. See plots below.

    Code:
     (* First call to get u01[t] *)
    result := 
     NDSolve[{Derivative[1][u][t] == (k4 - u[t])*k1 + (w[t] - u[t])*k2, 
       Derivative[1][w][t] == (u[t] - w[t])*k3, u[0] == k5, 
       w[0] == k5}, {u[t], w[t]}, {t, 0, k6}]; 
    u01[t_] = First[Evaluate[u[t] /. result]];
    (* now cycle functions into NDSolve *)
    myCurrentFunction[t] = u01[t];
    myFunctionTable = Table[
      result := 
       NDSolve[{Derivative[1][u][
           t] == (myCurrentFunction[t] - u[t])*k1 + (w[t] - u[t])*k2, 
         Derivative[1][w][t] == (u[t] - w[t])*k3, u[0] == k5, 
         w[0] == k5}, {u[t], w[t]}, {t, 0, k6}]; 
      myCurrentFunction[t] = First[Evaluate[u[t] /. result]];
      myCurrentFunction[t], {n, 1, 8}]
    myFunctionTable = Prepend[myFunctionTable, u01[t]]
    tablelength = Length[myFunctionTable]
    ms = Table[Plot[myFunctionTable[[i]], {t, 0, k6}], {i, 1, tablelength}]
    GraphicsArray[{{ms[[1]], ms[[2]], ms[[3]]}, {ms[[4]], ms[[5]], 
       ms[[6]]}, {ms[[7]], ms[[8]], ms[[9]]}}]
    Attached Thumbnails Attached Thumbnails Mathematica - ODE system solving problem-garray.jpg  
    Last edited by shawsend; November 26th 2009 at 10:16 AM.
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Newbie
    Joined
    Nov 2009
    Posts
    8
    Great! It works! Thank You Very Much, Shausend!
    Here's a piece of code providing the result which i have been trying to obtain for the last three days /ODE system was amplified with one more equation, but it doesn't influence on behavior of u[t] at all/:
    Code:
    k1=25;k2=0.3;k3=0.04;k4=0;k5=0.27;
    k6=0.98;k7=0.002;k8=0.03;k9=0.004;k10=100;nt=10
    result := NDSolve[{u'[t] == (k6 - u[
        t])*k1 + (v[t] - u[
            t])*k2, v'[t] == (u[t] - v[t])*k3, w'[t] == (k8 - w[t])*k1 + (u[t]*(
          1 - v[t]) + v[t]*(1 - u[t]) - w[t])*k4 + (2*v[t](1 - v[t]) - w[t])*k5, 
            u[0] == k7, v[0] == k7, w[0] == k9}, {u[t], v[t], w[t]}, {t, 0, k10}]
    uf[t_] = First[Evaluate[u[t] /. result]]
    wf[t_] = First[Evaluate[w[t] /. result]]
    pfr = Table[
        result := NDSolve[{u'[t] == (uf[t] - u[t])*k1 + (v[t] - u[t])*k2, v'[
          t] == (u[t] - v[t])*k3, w'[t] == (wf[t] - w[t])*k1 + (u[t]*(1 - v[
        t]) + v[t]*(1 - u[t]) - w[t])*k4 + (2*v[t](1 -
       v[t]) - w[t])*k5, u[0] == k7, v[0] == k7, w[0] == k9}, {u[t], v[t], w[
        t]}, {t, 0,
       k10}]; uf[t_] = First[Evaluate[
          u[t] /. result]]; wf[t_] = First[Evaluate[w[t] /. result]];
        {uf[t], wf[t]}, {n, 2, nt}]
    Plot[{0, uf[t], wf[t], 1}, {t, 0, k10}]
    output := Table[tr := n; {t, uf[tr], wf[tr]}, {n, 1, k10}]
     Export["output.dat", output]
    See the picture plotting in Mathematica. These curves are exactly the things I wanted to obtain.
    By the way, i tried to apply this algorithm to the "For" procedure, but it didn't work. I didn't know that Table is more convenient for such purposes in Mathematica.

    Not sure why those 6 and more are flat lining though
    Because the range of arguments in plots is strangely small.

    One small problem is still leaving: in "out" log there are too much lines with "Interpolating function..." appearing as a result of '=' operation inside the table. it is not good if number of steps nt is 500, all this inscriptions are unnecessary for me. How to disable this option in "Out"? I want to see only the plot graph.
    Last edited by LeChat; November 27th 2009 at 06:31 AM.
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Super Member
    Joined
    Aug 2008
    Posts
    903
    The semicolon at the end of a command disables output to the screen although I believe that was implemented in ver 6. So for me, if I just add a semicolon at the end of the table command, the output of all the interpolation messages is disabled. If that doesn't work in ver 5.4, there may be another way to disable the output to the screen.
    Follow Math Help Forum on Facebook and Google+

  9. #9
    Newbie
    Joined
    Nov 2009
    Posts
    8

    Thumbs up

    Ok, all problems are over now!

    Of course, I might see tutorials instead of asking such simple things here and thus wasting your time, but I haven't any time for it, sorry.

    And again, TYVM, Shausend! I've donated money to you a little.
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Replies: 5
    Last Post: August 4th 2011, 03:36 PM
  2. Replies: 3
    Last Post: May 3rd 2011, 11:03 PM
  3. Solve system in Mathematica in terms of parameters
    Posted in the Math Software Forum
    Replies: 6
    Last Post: February 25th 2011, 08:24 AM
  4. Not solvable system of linear equations in Mathematica
    Posted in the Math Software Forum
    Replies: 8
    Last Post: November 5th 2010, 11:22 AM
  5. Mathematica: Solving equations
    Posted in the Math Software Forum
    Replies: 7
    Last Post: March 12th 2010, 02:26 PM

Search Tags


/mathhelpforum @mathhelpforum