Maxima comes with a few well-implemented numerical ODE solvers — rk() and rkf45() — and I’ve written in a previous post about my variable stepsize implementation of the 2nd order backward differentiaion formula for stiff systems BDF2a().
Plotting the resulting numerical solutions — which is a list of lists of the form
[[indvar, depvar1, depvar2,… ], [indvar, depvar1, depvar2,… ], … ]
with each of the sublists representing the state of the system at a discrete time — requires a little familiarity with list functions like makelist() and map().
Here are a few little helper routines to make that quicker:
wxtimeplot()
This is mainly a teaching tool for my ordinary differential equations course, taking the output of the solver and, assuming there aren’t more than three dependent variables, plotting all the dependent variables vs the independent variable on the same set of axes, with a representation of the step size.
wxphaseplot2d,wxphaseplot3d
I wrote about these in another post
sol_points()
This is a more general purpose function. Given the output of the numerical ode solver, this constructs a points() subcommand ready to plug into draw2d(). This allows for the full generality of draw2d: colors, axis labels, point points,
Especially useful for plotting results from more than one call to the solver, such as when we need to see the effect of changing a parameter value:
Below I’ve pasted the code for the plots, and also for the interesting differential equations from Larter et al in Chaos, V.9 n.3 (1999)
(tw:(cosh((Vi-V3)/(2*V4))^(-1)),
m:(1/2)*(1+tanh((Vi-V1)/V2)),
w:(1/2)*(1+tanh((Vi-V3)/V4)),
aexe:aexe1*(1+tanh(Vi-V5)/V6),
ainh:ainh1*(1+tanh((Zi-V7)/V6))
)$
(rhs1:-gCa*m*(Vi-1)-gK*Wi*(Vi-ViK)-gL*(Vi-VL)+I-ainh*Zi,
rhs2:(phi*(w-Wi)/tw),
rhs3:b*(c*I+aexe*Vi)
)$
(V1:-.01,V2:0.15,V3:0,V4:.3,V5:0,V6:.6,V7:0,
gCa:1.1,gK:2,gL:0.5,VL:-.5,ViK:-.783,phi:.7,tw:1,
b:0.0809,c:0.22,I:0.316,aexe1:0.6899,ainh1:0.695)$
load(rkf45);
sol:rkf45([''rhs1,''rhs2,''rhs3],[Vi,Wi,Zi],[.1,.1,.1],[t,0,400],
report=true,absolute_tolerance=1e-8)$
sol_points(numsol,nth,mth):=points(map(nth,numsol),map(mth,numsol));
/* wxtimeplot takes the output of sol: rk45 and plots
up to 3 dependent variables
plus the scaled integration step size vs time*/
wxtimeplot(sol):=block(
[t0,t,tt,dt,big,dtbig],
t0:map(first,sol),
t:part(t0,allbut(1)),
tt:part(t0,allbut(length(t0))),
dt:t-tt,
dtbig:lmax(dt),
big:lmax(map(second,abs(sol))),
if is(equal(length(part(sol,1)),3)) then(
big: max(big,lmax(map(third,abs(sol)))),
wxdraw2d(point_type=6, key="y1",
points(makelist([p[1],p[2]],p,sol)),
color=red, key="y2",
points(makelist([p[1],p[3]],p,sol)),
color=magenta, key="dt",point_size=.2,
points(t,big*dt/dtbig/4),xlabel="t",ylabel=""))
elseif is(equal(length(part(sol,1)),2)) then
wxdraw2d(point_type=6, key="y1",
points(makelist([p[1],p[2]],p,sol)),
color=magenta, key="dt",point_size=.2,
points(t,big*dt/dtbig/4),xlabel="t",ylabel="")
elseif is(equal(length(part(sol,1)),4)) then(
big: max(big,lmax(map(third,abs(sol)))),
big: max(big,lmax(map(fourth,abs(sol)))),
wxdraw2d(point_type=6, key="y1",
points(makelist([p[1],p[2]],p,sol)),
color=red, key="y2",
points(makelist([p[1],p[3]],p,sol)),
color=green, key="y3",
points(makelist([p[1],p[4]],p,sol)),
color=magenta, key="dt",point_size=.2,
points(t,big*dt/dtbig/4),xlabel="t",ylabel=""))
);
The MaximaList 