Numerical Solutions of 2D and 3D ODE systems and phase plots

I most often use the built-in solvers in MATLAB (or the freeware alternative FreeMat) for numerical solutions of ODE systems.  They are easy to call and MATLAB graphical output is very flexible and of high quality.  That said, if I’m already working in Maxima to study stability behavior and equilibrium solutions, it’s fun to stay in the same window for  a numerical solution.

Here’s an example:  a damped, driven hard spring oscillator, written in phase space variables

\dot{y}=z, \dot{z}=-y^3-y-z/10 + 3/10\cos(t)

It’s straight-forward enough to run the Maxima numerical solver rk on this equation.  Extracting the numerical values of the two dependent variables is a little bit of a drag, requiring a call to makelist.  To make a plot in phase space, here’s a one-line function that is so easy to call:

/* wxphaseplot2d takes the output of s: rk([rhs1,rhs2]) 
   and  plots in phase space */
wxphaseplot2d(s):=wxplot2d([discrete,makelist([p[2],p[3]],p,s)],[xlabel,""],[ylabel,""]);

Here’s how it works on the above  example with a cute phase-space solution curve:

phaseheart

And here’s another famous phase curve—the solution of the Brusselator:

brusselatorphase

 

And can I just point out how analytically solving for the equilibrium solution has a very similar calling sequence — so easy that there’s no reason not to do it every time you want a numerical solution, and vice versa!

 

brusselator_eqsoln

 

Update!  Here’s a 3D phase plotter

/* wxphasplot3d takes the output of s: rk([rhs1,rhs2,rhs3]) and plots in phase space */

wxphaseplot3d(s):=wxdraw3d(point_type=0,point_size=.3,points_joined=true,points(makelist([p[2],p[3],p[4]],p,s)));

A classical example of 3D phase space trajectory curves, the Lorenz strange attractor system:

lorenz3D

Advertisements

2 thoughts on “Numerical Solutions of 2D and 3D ODE systems and phase plots”

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s