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
)$

Sponsored Post Learn from the experts: Create a successful blog with our brand new courseThe WordPress.com Blog

Are you new to blogging, and do you want step-by-step guidance on how to publish and grow your blog? Learn more about our new Blogging for Beginners course and get 50% off through December 10th.

WordPress.com is excited to announce our newest offering: a course just for beginning bloggers where you’ll learn everything you need to know about blogging from the most trusted experts in the industry. We have helped millions of blogs get up and running, we know what works, and we want you to to know everything we know. This course provides all the fundamental skills and inspiration you need to get your blog started, an interactive community forum, and content updated annually.

C’est une pipe for Maxima

 

I really love the pipe operator $>$ introduced into R by the package magrittr.  Not least because of the reference to the piece pictured here — “The Treachery of Images” by the Belgian surrealist painter René Magritte.

I don’t know yet whether this will be a really useful Maxima feature, but it is easy enough to make a simple pipe operation using infix():

infix("%>%",1,1); 
"%>%" (a,b):=b(a);

pipe1

Notice the 2nd and 3rd arguments of infix() control the precedence of operations (the left and right hand binding powers lbp and rbp).  We want to set those values lower than lbp and rbp for multiplication * in 2*%pi  and for for addition + in 1+%i (which by default in Maxima have lbp=180 and rbp=180).  Here’s the maxima document that gives the details about operators and binding power.

To make the pipe idea slightly more powerful so as to work with functions like integrate() and diff() which require more than one argument, we could try something like this, which makes use of the matlab-like find() I wrote about recently.

infix("%>%",1,1);
"%>%" (a,b):=block([i,bb],
 if atom(b) then 
     b(a)
 else (
     i:find(b=%P),
     bb:append(firstn(b,i[1]-1),[a],rest(b,i[1])),
     apply(first(bb),rest(bb,1))
     )
);

pipe2.PNG

Note that we’ve introduced the convention that the piped expression is referred to as %P, and that we call functions using a lisp-like construction where diff(%P,x) becomes [diff,%P,x].

A MATLAB-like find() function for Maxima

In MATLAB, you can find the indices of array entries that match a user-supplied condition with the function find()

Here’s a similar function in Maxima.  This makes liberal use of the functions sublist_indices()lambda(),  and several functions from the stringproc package: sequal()eval_string(), and simplode().  Oh, and also the shameless hack of using ascii(32) as the space character 🙂

Here are some examples.  The code is included at the bottom.

find

find(exp):=block(
 oo:op(exp),
 if sequal(oo,"=") or sequal(oo,">") or sequal(oo,">=") or sequal(oo,"<") or sequal(oo,"<=") or sequal(oo,"#") then 
   sublist_indices(first(exp),lambda([x1],eval_string(simplode([x1,op(exp),last(exp)]))))
 else 
   if sequal(oo,"and") or sequal(oo,"or") then(
     e1:first(exp),
     e2:second(exp),
     sublist_indices(first(e1),lambda([x1],
       eval_string(simplode([x1,op(e1),last(e1),ascii(32),oo,ascii(32),x1,op(e2),last(e2)]))))
   )
)$

Two little list utilities for MATLAB-like array indexing

I grew up with MATLAB, where extracting a subset of a vector V was easy as feeding a vector of indices into a vector.  For example entries i thru j could be had with V(i:j) and an more complicated index scheme could be accomplished with a vector of indices i_index and then V(i_index).

In Maxima, first(), last(), firstn(), rest(), and most generally makelist() allows for all that and more but with a little more cumbersome calling protocols.  Here are one-liners that achieve something like the two canonical MATLAB examples above:

indexby

A maxima function to replicate MATLAB diff() and an efficiency comparison

I needed a function to compute a list containing the difference between adjacent entries in another list, as diff() does for an array in MATLAB.

I first wrote  something using makelist() and found it very slow for large lists.  I then rewrote the function using rest() and found it much faster:

differ

 

Surface Integrals, Triple Integrals, and the Divergence Theorem of Gauss in Maxima

In earlier posts, I describe the Package of Maxima functions MATH214 for use in my multivariable calculus class, with applications to Greens  Theorem and Stokes Theorem.

Here we show how the  surface integral function integrateSurf() and triple integration  function integrate3()  (together with the divergence function div() )work on a Gauss’s Theorem example:

Gauss2

We integrate the parabolic surface and the circular base surface separately, and show their sum is equal to the triple integral of the divergence.

Gauss1

The functions above are included in the MATH214 package, but I list them below as well:

integrateSurf(F,S,uu,aa,bb,vv,cc,dd):=block(
 [F2],
 F2:psubst([x=S[1],y=S[2],z=S[3]],F),
 integrate(integrate(trigsimp(F2.cross(diff(S,uu),diff(S,vv))),uu,aa,bb),vv,cc,dd));

integrate3(F,xx,aa,bb,yy,cc,dd,zz,ee,ff):=block(
 integrate(integrate(integrate(F,xx,aa,bb),yy,cc,dd),zz,ee,ff));

div(f,x,y,z):=diff(f[1],x)+diff(f[2],y)+diff(f[3],z)$

cross(_u,_v):=[_u[2]*_v[3]-_u[3]*_v[2],_u[3]*_v[1]-_u[1]*_v[3],_u[1]*_v[2]-_u[2]*_v[1]]$

 

Path Integrals in the Plane, Double Integrals, and Greens Theorem in Maxima

In an earlier post I described the Maxima package MATH214 for use in my multivariable calculus class.  I’ve posted examples with applications to Gauss’s  Theorem and Stokes Theorem.

Here we take the double integration routine integrate2() and the 2D path integral integratePathv2() for a spin with a Green’s Theorem example from Stewart’s Calculus Concepts and Contexts:

Greens2

Greens1

And of course polar coordinates are nice too:

Greens3.PNG

The two functions used above are included in the MATH214 package, but I list them below as well.

integratePathv2(H,r,t,a,b):=block(
[H2],
H2:psubst([x=r[1],y=r[2]],H),
 integrate( trigsimp(H2.diff(r,t)),t,a,b)
);

integrate2(F,xx,aa,bb,yy,cc,dd):=block(
 integrate(integrate(F,xx,aa,bb),yy,cc,dd));

 

 

Surface Integrals and Stokes Theorem in Maxima

 

In an earlier post I detailed the Maxima functions contained  the MATH214 package for use in my multivariable calculus class.  The package at that link has now been updated with some further integration utilities:  integrate2() and integrate3() for double and triple integrals and integrateSurf() for surface integrals of vector fields in 3D.  I’ve posted examples with applications to Green’s  Theorem and Gauss’s Theorem.

Here’s a test drive of the surface integration function using a Stokes theorem example I found on the web:

Verify Stokes theorem for the surface S described by the paraboloid z=16-x^2-y^2 for z>=0

and the vector field

F =3yi+4zj-6xk

First the path integral of the vector field around the circular boundary of the surface using integratePathv3() from the MATH214 package

Stokes1

And also the surface integral using integrateSurf().  Notice that in the order of integration we specify in integrateSurf(),  (first  y then x) the surface normal vector computed with cross() in that same variable order points inward—the negative orientation.  We reverse the direction with an extra negative inside the surface integral.

Stokes2

Although they are included in the MATH214 package, here are the  functions used above:

integrateSurf(F,S,uu,aa,bb,vv,cc,dd):=block(
 [F2],
 F2:psubst([x=S[1],y=S[2],z=S[3]],F),
 integrate(integrate(trigsimp(F2.cross(diff(S,uu),diff(S,vv))),uu,aa,bb),vv,cc,dd));

cross(_u,_v):=[_u[2]*_v[3]-_u[3]*_v[2],_u[3]*_v[1]-_u[1]*_v[3],_u[1]*_v[2]-_u[2]*_v[1]]$

curl(f,x,y,z):=[ diff(f[3],y)-diff(f[2],z),curl(f,x,y,z):=[ diff(f[3],y)-diff(f[2],z),                     diff(f[1],z)-diff(f[3],x),   diff(f[2],x)-diff(f[1],y) ]$
integratePathv3(H,r,t,a,b):=block(
[H3],
H3:psubst([x=r[1],y=r[2],z=r[3]],H),
 integrate( trigsimp(H3.diff(r,t)),t,a,b)
);

 

 

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:

laplace1

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:

laplace2

So, I’ve written an alternative inverse Laplace function laplaceInv() that fixes that problem:

laplace3

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

laplace6

laplace7

A second-order equation with two impulsive forces:

laplace4

laplace5

The new Laplace transform Maxima function can be downloaded here.

Popularity of CAS programming languages: Maxima, Maple, Mathematica

A few months ago, inspired by the PYPL PopularitY of Programming Language,  I compared Google trends data for the 3M of CAS software:  Maxima, Maple, and Mathematica based searches of the form “<language> tutorial”.  The result was that Maxima seems to be slowly increasing in popularity  with about 20% of the interest in the 3M.

Today I saw another popularity metric: The TIOBE index.  Using their methodology of Google trends data for the search string <language> + programming, I have these  results for the proportion of searches among the 3M:

3MProgrammingProportions