Mathematica - ODE system solving problem

Printable View

November 25th 2009, 05:27 AM
LeChat

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!
November 25th 2009, 06:30 AM
shawsend

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}]
November 25th 2009, 01:41 PM
LeChat

Ok! (Happy) 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?
November 26th 2009, 04:17 AM
shawsend

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}]
November 26th 2009, 08:32 AM
LeChat

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 (Crying)
November 26th 2009, 09:00 AM
shawsend

1 Attachment(s)

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]]}}]
November 27th 2009, 03:24 AM
LeChat

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.

Quote:

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.
November 27th 2009, 08:33 AM
shawsend

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.
November 29th 2009, 01:07 PM
LeChat

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.