Maxima Language Code Translations

Have you seen these sites?

Hyperpolyglot.org  and Rosettacode.org  provide a valuable help for users who need to solve a problem in a new language:  Line-by-line comparisons  with lots of other languages.

For Maxima users, Hyperpolyglot’s computer-algebra system section  has side-by-side tables of comparisons with Mathematica, Maple, Sage and Numpy.

Rosettacode is arranged by task, with user-contributed solutions to common tasks in lots of languages.  I was motivated to write this post while working on a new Maxima version of an ordinary differential equations course I’ve taught for years using MATLAB.

 Here is Rosettacode’s section on the Euler Method  — a method for numerical approximation to the solution of a first order ordinary differential equation.

For the record, this is the version I decided to teach in my course:

(    /* Euler Method for  initial value problem
          y'=xy ,  y(0)=0.5   */ 
x0:0,
y0:.5,
h:.25,
nsteps:10,
xeuler:makelist(0,i,1,nsteps+1,1),
yeuler:makelist(0,i,1,nsteps+1,1),          
xeuler[1]:x0,
yeuler[1]:y0,
for i:1 thru nsteps do (
     xeuler[i+1]:xeuler[i]+h,
     yeuler[i+1]:yeuler[i]+h*xeuler[i]*yeuler[i]  
     )                                        
);
Advertisements

Checking the order of accuracy for the built-in ODE solver rk

The Maxima documentation says the numerical ODE solver in the command rk() is  fourth order accurate.  I ran a little test to confirm that.

Here’s the html export of the Maxima session, and a screenshot of the punchline:

RKorder4

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

Root of a nonlinear equation, as the last step of an ODE solution

What better use for a computer algebra system than a problem whose solution you know intuitively, and for which the paper-and-pencil work feels too daunting to start?  Here’s a problem that exploits what we know about the solutions of damped driven oscillators, making use of the usual ODE capabilities in Maxima, together with a nonlinear solve with find-root at the end.