Phase Plane Analysis of Linear Differential Equations

Description of behavior from coefficients of linear expansion about fixed point

Write a set of function that predict the behavior of the system  ((d  (Overscript[y,  ->]) )/(d t))= (({{a, b}, {c, d}})) (Overscript[y,  ->])  for real matrices.

Note that the eigenvalues are either real or they are complex conjugates (i.e., the eigenvalues are the solutions to  the quadratic equation)

Here is a function to print out a warning if the eigenvalues are inconsistent with a real coefficient matrix:

EigenConsistency[eval1_,  eval2_]  :=  
      If[And[And[Im[eval1]  ≠ 0, Im[eval2]  ≠  0], 
          Conjugate[eval1]  ≠  eval2], 
        Print["Coefficients are not real"]]

Here is a function that prints out information about stability based on the eigenvalues:

EigenStability[eval1_,  eval2_]  :=  
      If[And[Re[eval1] < 0, Re[eval2]  < 0], 
        Print["Stable and Attractive"], 
        If[Or[Re[eval1] > 0, Re[eval2] > 0], Print["Unstable"], 
          Print["Stable Orbits about Fixed Point"]]]

Here is a function that prints out information about the type of fixed point:

EigenTrajectory[eval1_,  eval2_]  :=  
      If[Or[Im[eval1]  ≠  0 , Im[eval2]  ≠  0], 
        Print["Circulation"], 
If[(s1 = 
                Sign[Re[eval1]])  ≠  (s2 = 
                Sign[Re[eval2]]), Print["Saddle"], 
          Print["Node"]]]

Here is a function that collects our other functions

EigenDescription[eval1_,  eval2_]  :=  
      Module[{}, EigenConsistency[eval1, eval2]; 
        EigenStability[eval1, eval2]; 
        EigenTrajectory[eval1, eval2]]

EigenDescription[ -1  + I, -1  - I]

EigenConsistency[1 +  I, 1]

Here is a function that prints out information about the direction of the eigenvectors

EigenDirector[eval_,  {ex_,  ey_}]  :=  
      Module[{theta  =  N[180*ArcTan[ex, ey]/π]}, 
        If[eval > 0,  
          Print["Unstable (λ=", 
            eval, ") direction is θ = ",  
            theta]]; 
If[eval < 0,  
          Print["Stable (λ=", 
            eval, ") direction is θ = ",  theta]]]

Here is a function that uses the coefficients to collect and print the physical information about the fixed point

LinearDescription[a_,  b_ ,  c_ ,  d_]  :=  
      Module[{esys, eval1, eval2, evec1, evec2}, 
esys = 
          Eigensystem[{{a, b}, {c, d}}]; eval1 = Chop[esys[[1, 1]]]; 
        eval2 = Chop[esys[[1, 2]]]; evec1  =  Chop[esys[[2, 1]]]; 
        evec2 = Chop[esys[[2, 2]]]; 
EigenDescription[
          eval1, eval2]; 
Print["Eigenvalues = ",  
          eval1,  " and ",  eval2]; 
If[
          And[Im[eval1]  ==  0,  Im[eval1]  ==  0],  
          EigenDirector[eval1, evec1];  EigenDirector[eval2, evec2]]; 
        esys]

LinearDescription[1,  .1, 1, -3]

Visualization of Trajectories near Fixed Points in the Plane

Functions to find the solutions for input coefficients (first one for specified initial point; the second one for initial points picked randomly

LinSol[a_,  b_ ,  c_ ,  d_  ,  xo_ ,  
        yo_]  :=  
{x[t], p[t]} /. 
        Flatten[DSolve[{x'[t]  ==  
                a x[t]  +  b p[t], p'[t]  ==  
                c x[t]  +  d p[t],  x[0]  ==  xo,  
              p[0]  ==   yo}, {x[t], p[t]}, t]]

LinSol[a_,  b_ ,  c_ ,  
        d_ ]  :=  
{x[t], p[t]} /. 
        Flatten[DSolve[{x'[t]  ==  
                a x[t]  +  b p[t], p'[t]  ==  
                c x[t]  +  d p[t],  
              x[0]  ==  Random[Real, {-10, 10}],  
              p[0]  ==   Random[Real, {-10, 10}]}, {x[t], p[t]}, 
            t]]

LinSol[3, 2, -3, 2]

Example of a plot for a single directory

CritPointPlot[a_,  b_ ,  c_ ,  d_ ]  := 
      ParametricPlot[Evaluate[LinSol[a, b, c, d]], {t, 0, 20}, 
        PlotRange -> {{-15, 15}, {-15, 15}}]

CritPointPlot[-1, - .2, 1, -2]

This is not terribly informative.  Construct a function that takes input and plots many paths with arrows that describe the trajectory.  Print out physical interpretations as well.  Supply an argument with  a default definition that indicates the number of trajectories to plot.

<< Graphics`Arrow`

CritPointPlotPointsMany[a_,  b_ ,  c_ ,  
        d_ ,   howmany_:  100]  := 
      Module[{esys, eval1, eval2, funcs , data, lendata, gstack, lp, rstack, 
          lstack, magrange, rstackmag, lstackmag}, 
esys = 
          LinearDescription[a, b, c, d]; Print[esys]; 
        eval1 =  esys[[1, 1]];  
        eval2 =  esys[[1, 2]]; 
funcs = 
          Table[Chop[LinSol[a, b, c, d]], {i, howmany}];  lstack  = {}; 
        rstack = {}; 
        For[imany = 1, imany  ≤  howmany, imany++, 
          data = Chop[
              Table[Evaluate[funcs[[imany]] /. t -> itime], {itime, 
                  0, 20,  .1}]]; lendata  =  Length[data]; 
          gstack = {Hue[0], 
              Arrow[data[[1]], data[[2]], 
                HeadScaling -> Absolute]}; 
For[
            iend  =  4,  iend  ≤   lendata/2, 
            iend  +=  8,  AppendTo[gstack, Hue[iend*0.66*2/lendata]]; 
            AppendTo[gstack, 
              Arrow[data[[iend - 1]], data[[iend]], 
                HeadScaling -> Absolute]]]; 
          lp = ListPlot[data, PlotJoined -> True, AspectRatio -> 1, 
              PlotRange -> {{-15, 15}, {-15, 15}}, 
              PlotStyle -> Hue[Random[]], 
              DisplayFunction -> Identity]; 
          lstack = Show[lp, Graphics[gstack], lstack]; 
          rstack  =  Show[lp, rstack]]; 
If[
          Or[Re[eval1] > 0, Re[eval2] >  0], 
          magrange = {{-60, 60}, {-60, 60}}, 
          magrange =  {{-1, 1}, {-1, 
                1}}]; 
rstackmag  =  
          Show[rstack, 
            PlotRange ->  
              magrange]; 
lstackmag  =   Show[lstack, 
            PlotRange ->  magrange]; 
Show[
          GraphicsArray[{{lstack, rstack}, {lstackmag, rstackmag}}], 
          ImageSize -> 1000, 
          DisplayFunction -> $DisplayFunction]]

CritPointPlotPointsMany[-1, 0.25, 0.75, -1]

CritPointPlotPointsMany[1.25, -0.25, -0.75, 0.75, 
      200]

CritPointPlotPointsMany[1.25, 0.25, 0.75, -1.75,  200]

CritPointPlotPointsMany[-0.25, -2.0, 3, 1.5]

CritPointPlotPointsMany[-0.25, -2.0, 2, -1.5]

CritPointPlotPointsMany[-0.25, -2.0, 2, 0.25]