3
$\begingroup$

I want to plot a system of three stiff ODEs (Oregonator model). That model describes a chemical oscillator. I do not have much experience with plotting ODEs, but I have obtained a Mathematica file that does almost the same as I want, but only with a simplified version of the Oregonator. That simplified version has two ODEs, while my system has three ODEs.

I have the following code:

ε=4*10^-2;
δ=4*10^-4;
q=8*10^-4;
f=1;

{xsol,zsol}=NDSolve[{εx'[t]==x[t](1-x[t])+(f(q-x[t])z[t])/(q+x[t]),z'[t]==x[t]-z[t],
                     x[0]==.00012,z[0]==.00576},{x,z},{t,0,40},MaxSteps->Infinity]

This system describes the following model:

$$ \epsilon \frac{dx}{dt}=x(1-x)+\frac{f(q-z)}{q+x}z $$ $$ \delta \frac{dz}{dt}=x-z $$

How do I modify the above code, so that it solves the below system?

ODE_SYSTEM

I plot the solutions like this:

Plot[Evaluate[x[t]/.xsol],{t,0,40},PlotRange->All,PlotStyle->{Thick,Blue}]

I already tried adding ysol and y[t] and y[0], but that didn't work for me. Probably I did something wrong, but I'm really stuck. The initial value is y[0]=0.375.

How should I modify the code with the NDSolve? Thanks a lot!

REQUEST: it is appreciated if sb replaces the images of the system with the the correct MathJax-notation. I don't have experience with MathJax.

$\endgroup$
  • $\begingroup$ Should the last term in the first equation of the 3-species system be $x(1-x)$ instead of $x(1-z)$? $\endgroup$ – Michael E2 May 23 at 15:45
  • $\begingroup$ @MichaelE2 That's indeed true. Sharp eye do you have ;). Thanks for being so attentive leaving a comment! I'll change the question and the answer so that it represents the correct system. $\endgroup$ – ralphsmit May 25 at 7:03
7
$\begingroup$

How do I modify the above code, so that it solves the below system?

You did not have IC for $y(t)$, so I made one up.

Clear["Global`*"];
ε = 4*10^-2;
δ = 4*10^-4;
q = 8*10^-4;
f = 1;

ode1 = ε x'[t] == q y[t] - x[t]*y[t] + x[t] (1 - x[t]);
ode2 = δ y'[t] == -q y[t] - x[t]* y[t] + f*z[t];
ode3 = z'[t] == x[t] - z[t];

ic = {x[0] == .00012, y[0] == .375, z[0] == .00576};

{xSol, ySol, zSol} = NDSolveValue[{ode1, ode2, ode3, ic}, {x, y, z}, {t, 0, 40}]

To plot, one way

p[f_] := Plot[f, {t, 0, 40}, PlotRange -> All];
Grid[{{p[xSol[t]], p[ySol[t]], p[zSol[t]]}}, Frame -> All]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Hi @Nasser thanks for also updating the images and of course the answer in the first place. In the first image for x[t], some parts of the graph are not displayed. Do you know how that comes and how I can display the first graph completely? So that I have an idea what the maximum is for x[t]. $\endgroup$ – ralphsmit May 25 at 9:54
  • 2
    $\begingroup$ @ralphsmit updated. You can use PlotRange->All without this option, Mathematica makes its own mind what is a "good" range. $\endgroup$ – Nasser May 25 at 9:56
  • $\begingroup$ @Nassar Awesome! Thanks for telling met that! $\endgroup$ – ralphsmit May 25 at 9:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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