Nonlinear Chemical Dynamics
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

Governing equation can be found in:

G. Nicolis, Dynamique Chimique, Thermodynamique, cinétique et mécanique statistique, Dunod, Paris 2005. 

The author would like to acknowledge the help of Professor Brian Higgins, Chemical Engineering Department, U C Davis who kindly provided notebooks on continuation methods.

We find the steady states for a fixed value of λ=1.  We get the famous S-shape curve. Then, we look for the steady states for a fixed value of μ=0.2. We get two disjoint braches (one branch has a turning point and gives two steady states). Here, we have used the built-in function of Mathematica called: ImplicitPlot.

λ = 1; μ =.

Needs["Graphics`ImplicitPlot`"]

plt1 = 
      ImplicitPlot[-z^3 + λ z + μ == 
          0 , {μ, -10, 10}, {z, -5, 5}, AspectRatio ->  .7, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotPoints ->  1000]

λ =. ; μ = 0.2;

plt2 = 
      ImplicitPlot[-z^3 + λ z + μ == 
          0 , {λ, -10, 10}, {z, -5, 5}, 
        AspectRatio ->  .7, PlotStyle -> RGBColor[1, 0, 0], 
        PlotPoints ->  1000]

In this section, we use arc length continuation methods.

λ = 1; μ =. ;

eq1 = -z^3 + λ z + μ == 0 /. 
          z -> z[s] /. μ ->  μ[s]

eq2=μ'[s]^2== 1-((z'[s])^2)

init1 = μ[0] == 0

Solve[-z^3 + λ z + μ == 0 /. μ ->  0, 
      z]

init2 = z[0] == -1

sol=NDSolve[{eq1,eq2,init1,init2,μ'[0]== 1},{z[s], μ[s]},{s, -100, 100},SolveDelayed -> False]

plt3 = 
        ParametricPlot[Evaluate[{μ[s], z[s]} /. sol], {s, -10, 10},  
          PlotStyle -> RGBColor[0, 0, 1], 
          PlotRange -> All];

Show[plt1, plt3]

λ =. ; μ = 0.2;

eq1 = 
      D[-z^3 + λ z + μ == 0 /. 
            z -> z[s] /. λ ->  λ[s], 
        s]

eq2=λ'[s]^2== 1-((z'[s])^2)

init1 = λ[0] == 0

Solve[-z^3 + λ z + μ == 0 /. λ ->  
          0, z]

init2 = z[0] == 0.584804

sol = 
      NDSolve[{eq1, eq2, init1, init2}, {z[s], λ[s]}, {s, -100, 
          100}, SolveDelayed -> False]

plt4 = 
        ParametricPlot[
          Evaluate[{λ[s], z[s]} /. sol], {s, -10, 10},  
          PlotStyle -> RGBColor[0, 0, 1], 
          PlotRange -> All];

Show[plt2, plt4]

init1 = λ[0] == 1

Solve[-z^3 + λ z + μ == 0 /. λ ->  
          1, z]

init2 = z[0] == -0.878885

sol = 
      NDSolve[{eq1, eq2, init1, init2}, {z[s], λ[s]}, {s, -100, 
          100}, SolveDelayed -> False]

plt5 = 
        ParametricPlot[
          Evaluate[{λ[s], z[s]} /. sol], {s, -10, 10},  
          PlotStyle -> RGBColor[0, 0, 1], 
          PlotRange -> All];

Show[plt2, plt5]