## Multiple Plots in an Animated GIF in Maxima

Here it is: Roughly the animation I wanted to create. A harmonic oscillator, with the corresponding trajectory in phase space and a bead-on-wire projection of the trajectory onto the potential energy curve.

I wanted to insert this in a powerpoint slide I would be using for an online lecture, so an animated GIF seemed like a good option. I started with this example from Dirk Mittler’s blog.

What I really wanted to do was include several axes like in these examples from my earlier post at TheMaximalist, which created multiple scenes using several calls to gr2d() and then drew the scenes in a single call to draw(). The trouble I encountered was that the collection time-synchronized multiple scenes weren’t rendered simultaneously in the animated GIF — the animation flickered back and forth between the scenes.

So I turned off the axes and made my own triplet of axes in a single scene using parametric lines and including the needed vertical offsets.

One key thing that took me a while to find: In my Windows installation (Maxima 5.44 and wxMaxima 20.06.6), no matter what I specified for a path in the file_name argument to draw(), the animated GIF file ended up getting written in C:\maxima-5.44.0

I’ve pasted the code below.

scenes: []$for i:0 thru 61 do ( scenes: append(scenes, [gr2d(proportional_axes=xy, noframe, ytics=false,xtics=false,nticks=100, /* oscillator 2 units below */ color=black, parametric(t,-2,t,-1.2,1.2), /* phase space 2 units above */ color=black, parametric(t,0,t,-1.2,1.2), parametric(0,t,t,-1.2,1.2), color=blue, parametric(cos(t),sin(t),t,0,2*%pi), /* potential */ color=black, parametric(t,2,t,-1.2,1.2), parametric(0,2+t,t,-.1,1), color=blue, explicit(2+x^2,x,-1.2,1.2), color=red,point_type=filled_circle, points([cos(i/10),cos(i/10),cos(i/10)],[-2,-sin(i/10),2+cos(i/10)^2]))] ) )$

draw(
delay = 10,
file_name = "oscillator_phasespace_potential",
terminal = 'animated_gif,
scenes
)$ ## An improved Maxima function for inverse Laplace transform Maxima has a fairly serviceable Laplace transform utility built-in. Here’s an example from the popular ordinary differential equations book by Blanchard, Devaney and Hall: Trouble arises when we look at discontinuous forcing functions, which is especially sad because it seems to me that’s what makes it worthwhile to spend time with Laplace transforms. In particular, the inverse transform function ilt() fails on the Heaviside Function and Dirac Delta, even though the built-in laplace() treats them correctly: So, I’ve written an alternative inverse Laplace function laplaceInv() that fixes that problem: Here are a few differential equation solutions to show how the new function behaves: A second-order linear equation with a constant forcing function that vanishes at $t=7$ A second-order equation with two impulsive forces: The new Laplace transform Maxima function can be downloaded here. ## logabs, logarc: How to make integrate() return what you expect ## logarc In the back of my calculus book there is a table of famous integrals. Here’s integral number 21 in that table: From Maxima integrate(), I get What’s going on? Both forms give a workable antiderivative for the original integrand: Furthermore, we believe that both forms are correct because of this helpful identity for hyperbolic sine: ${\rm asinh}(z)=\ln(z+\sqrt{1+z^2}).$ Turns out (thanks to a Barton Willis for pointing me in the right direction) there’s a variable logarc that we can set to make Maxima return the logarithmic form instead of hyperbolic sine: I haven’t yet encountered cases where this would be a bad idea in general, but I’ll update this if I do. ## logabs In the first week of my differential equations course, we study methods of direct integration and separation of variables. I like to emphasize that the absolute values can lend an extra degree of generality to solutions with antiderivatives of the form $\int \frac{1}{u}\;du = \ln |u|+ C.$ As an example, for the initial value problem $y' = \frac{x}{1-x^2}$, $y(0)=1$, it is convenient for treating all possible initial conditions ($x \ne \pm 1$) in one step to use the antiderivative $y = -\frac{1}{2} \ln | 1-x^2 | + C.$ However, Maxima omits the absolute values. For this case, we could consider only the needed interval $-1, but still… Turns out we can set the Maxima variable logabs to make integrate() include absolute values in cases like this: But then later in the course, I saw that logabs also impacts the Ordinary Differential Equation solver ode2(). I encountered an example for which Maxima, in particular solve() applied to expressions involving absolute value, didn’t do what I wanted with logabs:true For the logistic equation $\frac{dP}{dt} = kP\left ( 1-\frac{P}{P_c} \right )$, $P(0)=P_0$ we expect that by separating variables we can obtain the solution $P(t) = \frac{P_0 P_c}{(P_c-P_0)e^{-kt}+ P_0}.$ Here’s what happens when we use ode2() with and without logabs:true: ## Student Projects in Maxima Vol 1, Num 4: ODE Systems for Epidemiological Models Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences. Volume 1, Number 4 (2017) is devoted to Systems of Ordinary Differential Equations for SIR models in epidemiology. Contents: Emenheiser, Anna, Virus Dynamics and Drug Therapy Radermacher, Erin, A Model of the 2014 Ebola Outbreak Rafai, Sagar, Analysis of Communicable Disease Models ## Student Projects in Maxima Vol 1, Num 3: ODE Systems for Interacting Population Models Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences. Volume 1, Number 3 (2017) is devoted to Systems of Ordinary Differential Equations for non-chaotic predator-prey and other interacting population models Contents: Bhullar, Abhjeet, Lions, Wildebeest and Zebras Goodin Aunic, Parasitic Nematodes in Grouse DeFore, Brandon, Breeding Suppression in Predator-Prey Jerez, Emilio, Predator-Prey with Hawk and Dove Tactics Bontrager, Eric, Predator-Prey with Prey Switching Beatty, Ethan, Analysis of Logistic Growth Models Rice, Gabriel, Pharmacokinetic Mechanism of Ethanol-Benzodiazepine Interactions Kay, Ian, Radioactive Isotopes Wile, Jessica, Ebola in Western Lowlands Gorillas Bailey, John, Kinetic Modeling for Interconnected Reactions Piet, Joe, Elephant and Tree Population Dynamics Kim, Judy, A Model for West Nile Virus Kamalaldin, Kamalaldin, Infected Prey in Polluted Environment Park, Kayla, Ratio-Dependent Predator Prey Model Lundy, Liam, Predator Prey in a Single Species with Cannibalism Orwin, Michael, Interactions of Model Neurons Schultz, Pete, Photosynthetic Oscillations Wadhwa, Raoul, Dynamics in an Inflammation Model Del Olmo, Ricardo, Population Growth and Technological Change Network Model of Cocaine Traffic in Spain Edited Kill, Sean, A Logistic Model with Varying Carrying Capacity McFadden-Keesling, Sophia, Predator Prey Model with a Refuge Samson, Tanush, Lions and Wildebeest Model Morales, Zach, Game Theory Models ## Student Projects in Maxima Vol 1, Num 1: Systems of ODEs and Chaos Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences. Volume 1, Number 1 (2017) is devoted to Chaotic Systems of Ordinary Differential Equations. Contents: Thornburg, Eric, Simple Pendulum and Chaos Bhimani, Kevin, Chaos in a 3-Species Food Chain York, Lily, A Lattice Model of Epilepsy Andrews, Steven, 3D Chaotic Model Rutledge, Tim, A Model of Neuronal Bursting ## Student Projects in Maxima Vol 1, Num 2: Systems of ODEs with Impulses and Switching Functions Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences. Volume 1, Number 2 (2017) is devoted to Systems of Ordinary Differential Equations with time-dependent forcing terms, non-continuous inputs, and forcing terms that require knowledge of a past state of the system. Contents: Chumley, Qynce, An Office Heating and Cooling Model Williams, Nick, Baseball Pitch Dynamics Rizzolo, Skylar, BAC Model for alcohol consumption Barth, Eric, A model of bladder bacteria proliferation in prostate disease Barth, Eric, A Model of Pulse Vaccination Strategyi ## A Package of Maxima Utilities for my Ordinary Differential Equations Course: MATH280.mac I’ve put together a collection of functions — some direct quotes of other contributed functions, some renamed or repackaged, and some newly implemented — for various needed tasks in my undergraduate ordinary differential equations course. I’ve written elsewhere about the Backward Difference Formula implementations, the phase space visualization functions, the matrix extractors, and the numerical solutions plotters. The package includes my home-grown help utility. You can download the package MATH280.mac …and if you’re interested, here’s my multivariable calculus package MATH214.mac MATH280.mac contains: wxphaseplot2d(s) wxphaseplot3d(s) phaseplot3d(s) wxtimeplot(s) plotdf(rhs) wxdrawdf(rhs) sol_points(numsol,nth,mth) rkf45(oderhs,yvar,y0,t_interval) BDF2(oderhs,yvar,y0,t_interval) BDF2a(oderhs,yvar,y0,t_interval) odesolve(eqn,depvar,indvar) ic1(sol,xeqn,yeqn) ic2(sol,xeqn,yeqn,dyeqn) eigU(z) eigdiag(z) clear() - - for any of the above functions, help(function_name) returns help lines for function_name - Last Modified 5:00 PM 3/27/2017  ## Plotting Numerical Solutions of Ordinary Differential Equations in Maxima 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)$

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=""))
);