Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote -3 down vote favorite

I'm trying to solve this set of non linear differential equations, but apparently it won't solve. I really tried all the tweaks I could think of. Still no progress... So now I turn to you, can you help me? the code below should be ready to put into Mathematica. These equations describe the Belousov-Zhabotinsky reaction in chemistry.. 

(*Define the differential equations*)

deq1 := x'[t] == p[1] A y[t] - p[2] x[t] y[t] + p[3] B x[t] - 2 p[4] x[t]^2;

deq2 := y'[t] == -p[1] A y[t] - p[2] x[t] y[t] + f p[5] z[t] ;

deq3 := z'[t] == 2 p[3] B x[t] - p[5] z[t];


(*Replace A, B and all the p values, and set up the equations as a list*)

deqsys = {deq1, deq2, deq3 } /. {p[1] -> 1.7, p[2] -> 1.8*10^9, 
p[3] -> 9.0*10^3, p[4] -> 4*10^7, p[5] -> 0.5, f -> 2, B -> 0.6, 
A -> 0.6};

(*Solve the set of equations, and save them as the variable Sol*) 

Sol = NDSolve[deqsys, {x[t], y[t], z[t]}, t]
share|improve this question
1  
It's a non-linear system. There's no reason to think that it can be solved in closed form, except in fairly special cases. You need to look into NDSolve. There are loads of NDSolve examples in the documentation and on this site. Here's one: mathematica.stackexchange.com/questions/13072/… – Mark McClure Dec 14 '12 at 13:18
@MarkMcClure yes, I know it's non-linear, but I also tried the NDSolve, which still got me nowhere... – TehHO Dec 14 '12 at 13:25
@TehHO Then I suggest you rephrase your question in terms of NDSolve, since we don't expect DSolve to work anyway. – Mark McClure Dec 14 '12 at 13:27
@MarkMcClure, I just added an n. Doing this gives me an error NDSolve::ndnco: "The number of constraints (0) (initial conditions) is not equal to the total differential order of the system plus the number of discrete variables (3)." When i add initial conditions I get a really odd result, something with interpolating functions..? (not copyable for some reason..) – TehHO Dec 14 '12 at 13:31
2  
That's the solution – ssch Dec 14 '12 at 13:46
show 2 more comments

closed as too localized by acl, tkott, Sjoerd C. de Vries, Artes, rm -rf♦ Dec 15 '12 at 14:53

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, see the FAQ.

1 Answer

up vote 5 down vote

Why not something like this (just an example) : 

deq1 := x'[t] == 
  p[1] A y[t] - p[2] x[t] y[t] + p[3] B x[t] - 2 p[4] x[t]^2;
deq2 := y'[t] == -p[1] A y[t] - p[2] x[t] y[t] + f p[5] z[t];
deq3 := z'[t] == 2 p[3] B x[t] - p[5] z[t];
deqsys = {x[0] == 1.0, y[0] == 1.2, z[0] == 1.3, deq1, deq2, deq3} /. 
  {p[1] -> 1.7, p[2] -> 1.8*10^9, p[3] -> 9.0*10^3, 
   p[4] -> 4*10^7, p[5] -> 0.5, f -> 2, B -> 0.6, A -> 0.6} ;
Sol = NDSolve[deqsys, {x[t], y[t], z[t]}, {t, 0, 1}]

Plot[x[t] /. Sol , {t, 0, 1}]

enter image description here

share|improve this answer
Nice - and thanks for the contribution! I formatted your code by adding spaces at the beginnings of the lines and added a picture. If you examine the scale of the picture, you'll notice that the x range is essentially constant, at least with these initial inputs. You might try ParametricPlot3D to visualize the whole thing. – Mark McClure Dec 14 '12 at 15:32
like ParametricPlot[{y[t], z[t]} /. Sol // Evaluate, {t, 0.001, 1}] ? – chris Dec 15 '12 at 13:56

Not the answer you're looking for? Browse other questions tagged or ask your own question.