Complex roots of unity and the map() function

Here’s a little application from Complex Variables class that uses map() to repeatedly apply Maxima functions to each element in a list.

We want to find the complex 6th roots of 1 — there should be six such complex numbers z that solve the equation

  z^6=1

sixthrootseqn

Notice the solutions are returned as equations…this is sometimes really handy, but for now we want to isolate the right hand side of each equation. The command rhs(), when applied to an equation, returns just the right side. We can apply the rhs command to all the equations contained in our solution using the Maxima command map():

sixthrootsmap

Finally, with the real and imaginary parts extracted, we can plot the six roots

sixthrootsplot

 

An odd gnuplot error and a work-around using maxima_tempdir

***Update***********************

I also have this problem in the latest windows installation—on the same machine that long filenames were never a problem before in Maxima.  To fix it choose the sbcl version of maxima, as described at My Windows Installation.

*********************************

I have a student running wxMaxima on his windows laptop.  He has partitioned his hard drive so he can run Windows 10 as well as other operating systems.

When he attempts any plotting command, he gets an error:

gnuplot_error

The part of the error message with directory name TANUSH~1 says to me that there is some problem with the long directory name, possibly also a problem caused by a space in that directory name, possibly due to the way those new partitions are formatted.

As a workaround, we set the maxima option maxima_tempdir to point to another directory and all worked as expected:

maxima_tempdir

Of course, that new directory will now fill up with files with names like maxout_xxxx.dat and maxout_xxxx.gnuplot, but that was already happening in the default directory.  A natural feature for wxMaxima to include in future releases would be to clean up the directory by deleting all those files when the program is stopped…

Speeding Up Maxima Functions with compile()

In my differential equations class, we implement the first-order-accurate Euler method for numerical approximation of a solution of the initial value problem y' = f(x,y) , y(x_0)=y_0.

eulerloop

We implement this in an appealing way — from an algorithm translation point of view — but the inner loop with all that list indexing takes a lot of time.   To verify the accuracy of the method, we run the function with more and more successively smaller steps.

eulermethod(f,x0,y0,h,nsteps):=block(
 [i,x,y], 
 x:makelist(0,i,1,nsteps+1,1),
 y:makelist(0,i,1,nsteps+1,1),
 x[1]:x0,
 y[1]:y0,
 for i:1 thru nsteps do (
    x[i+1]:x[i]+h,
    y[i+1]:y[i]+h*ev(f,x=x[i],y=y[i])
 ),
 return ([x,y])
)$

Below I’ve used the time() function so you can see how simply running compile() on the function makes it run lots faster for ten thousand steps!

compileeuler

The Maxima documentation says this:

Function: compile
compile (f_1, …, f_n)
compile (functions)
compile (all)

Translates Maxima functions f_1, …, f_n into Lisp, evaluates the Lisp translations, and calls the Lisp function COMPILE on each translated function. compile returns a list of the names of the compiled functions.

compile (all) or compile (functions) compiles all user-defined functions.

compile quotes its arguments; the quote-quote operator '' defeats quotation.

Compiling a function to native code can mean a big increase in speed and might cause the memory footprint to reduce drastically. Code tends to be especially effective when the flexibility it needs to provide is limited. If compilation doesn’t provide the speed that is needed a few ways to limit the code’s functionality are the following:

  • If the function accesses global variables the complexity of the function can be drastically be reduced by limiting these variables to one data type, for example using mode_declare or a statement like the following one: put(x_1, bigfloat, numerical_type)
  • The compiler might warn about undeclared variables if text could either be a named option to a command or (if they are assigned a value to) the name of a variable. Prepending the option with a single quote ' tells the compiler that the text is meant as an option.

 

Here are the specs of the machine I ran these timings on:

processorinfo

 

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]  
     )                                        
);

A home-grown help utility for my Maxima functions

I have written a number of Maxima functions for my students.   I really miss the MATLAB convention by which help(functionname) prints to the console the leading comment lines in the function.  That is really easy for the developer, because the documentation lines within the function also serve as the user documentation in the console.

Maxima’s built-in documentation capabilities ?, ??, example(), demo(), and describe() are, as far as I can tell, unwieldy at best and unsuited for portable documentation to accompany user-written functions.

So, I’ve implemented my own home-grown help utility that can travel along with my *.mac files.  The scheme is:

  • Include in each *.mac file a simple help() utility that when invoked on the maxima command line as help(functionname) calls the ancillary documentation utility functionname_help() contained in the *mac file.  I’ll include this same little utility in all my *.mac files for consistency.
  • This requires that in addition to the functions themselves, the author of the functions include for each function in the *.mac file a help utility with the name convention functionname_help().
  • Upon loading of the *.mac file, I display an invitation to see the contents of the *.mac file using help(macFileName), which again requires an ancillary macFileName_help() function.  The resulting contents message also instructs in the use of the help() utility for functions.

Not as easy to implement as the MATLAB convention, but very workable from the point of view of my students as they learn to use these custom functions.

Below is an example for the ComplexVariables.mac, (which I’ll post in the next few days.)  I plan to use this as a template in all my *.mac files and I really hope the practice catches on in the Maxima community!

/*Complex Variables.mac 
Contents for use with help() */
ComplexVariables_help():=block(
printf(true,
"ComplexVariables.mac contains:
 factorC(f,z)
 partfracC(f,z)
 residueC(f,z,a)
 integratePathC(f,r,t,a,b)
 limitC(f,z,a)
 absC(exp)
 -
 -
 for any of the above functions,
 help(function_name) returns help lines for function_name"
)
);
/* The help utility */
help(_input):=block(
eval_string(sconcat(string(_input),"_help()")
)
);
/* the following is displayed upon loading */
disp("ComplexVariables.mac: help(ComplexVariables) for contents");

3D Printing and Maxima

In response to several nearly simultaneous queries, I’ve been working this week on generating a 3D  plot  and creating a mesh file from within Maxima  that is suitable for exporting to a 3D printer.

As a guide, I used this 2012 article from the Mathematical Intelligencer by Henry Segerman.  The idea is to draw a surface using traces made from parametric tubes, then export that to a .PLY format file.

Here’s my prototype for the process of generating the surface, reading in the resulting gnuplot data file, identifying faces appropriately, and writing into a .PLY file.

That isn’t yet ready for the 3D printer, but I could load it into the free program MeshLab, and from there convert to .STL

Here’s my figure from Maxima:

pringle

And the resulting .PLY file loaded into MeshLab

 

 

saddle_meshdraw

Here’s a few seconds of video from the printing process

And 2.5 hours later, the low resolution first-try 3D printed object:

3dprintverision01

Visualizing Complex Functions by their images over concentric circles

A few weeks ago I posted my Maxima code for visualizing complex function by their images over a grid of line segments in the complex plane.  Here are a few images from an experiment of plotting the image of complex function over concentric circles, centered at the origin, traversed in the counterclockwise direction.  The first image below shows the identity map.  The colors are in ROYGBIV rainbow order, as noted in the legend, which gives the radius of the pre-image circles.  The second image is for the function f(z)=(z-1/2)(z-1).  Notice for the smallest radius pre-image, the image under f is largely unchanged.  As radius grows the image forms a loop then doubles on itself.  Here’s a fun insight:  that looping behavior occurs at radius .75, a circle containing a point z at which f'(z)=0.  More on this later, but its exciting to find a geometric significance of a zero of the derivative, considering how much importance such points have in real calculus.  Finally the third figure shows really interesting behavior for the function f(z)=\frac{1}{(z-1/2)(z-1)}.  Notice again the looping brought about by points at which the derivative is zero.  Notice also the change in direction of the image each time we cross over a pole—that is a root of the denominator.  The maxima code is given below.

/* wxdrawcircleC
 plot the traces of complex valued function _f
 along concentric circles |z|=r for 0<r <=1
 blue - images of _f along horizontal lines
 red - images of _f along vertical lines */
wxdrawcircleC(_f):=block(
[fx,fy,fxt,fyt,pl,j,ncirc,r,theta_shift,lm,l,dx,dy,dxy,v,nt],
ratprint:false,
fx:realpart(rectform(subst(z=x+%i*y,_f))),
fy:imagpart(rectform(subst(z=x+%i*y,_f))),
pl:[],
nt:200,
r:2,
ncirc:7,
theta_shift:%pi/20,
/* get the max modulus for computing length of arrows */
lm:0,
for n:1 thru ncirc do (
l:cons(lm,makelist(cabs(subst(z=r*n/7*exp(%i*t),_f)),t,0,2*%pi,.5)),
 lm:lmax(l)
 ),
 for j:1 thru ncirc do(
 /*first the images of the circles */
 fxt:subst([x=r*cos(t)*j/ncirc,y=r*sin(t)*j/ncirc],fx),
 fyt:subst([x=r*cos(t)*j/ncirc,y=r*sin(t)*j/ncirc],fy),
 pl:cons([parametric(fxt,fyt,t,0,2*%pi)] , pl),
 pl:cons([color=colorlist(j)],pl),
 pl:cons([key=string(r*j/ncirc)],pl),
 /* now the arrows */
 dx:subst(t=3/nt,fxt)-subst(t=0,fxt),
 dy:subst(t=3/nt,fyt)-subst(t=0,fyt),
 dxy:sqrt(dx^2+dy^2),
 v:vector([subst(t=0,fxt),subst(t=0,fyt)],[lm*dx/dxy/20,lm*dy/dxy/20]),
 pl:cons([v],pl),
 pl:cons([head_length=lm,head_angle=1],pl),
 pl:cons([color=colorlist(j)],pl),
 pl:cons([key=""],pl)
 ),
pl:cons([xlabel="",ylabel=""],pl),
pl:cons([nticks=nt],pl),
pl:cons([grid=on],pl),
apply(wxdraw2d,pl)
);

Visualizing Complex Functions with Conformal Mapping

For the thoroughly modern calculus student, an introduction to Complex Variables is all the more daunting because we don’t have the kind of geometric intuition-building machinery available to us for functions f: C \rightarrow C as we did for real-valued functions in calculus.  Not that attempts haven’t been made….  In the coming weeks, I want to work on several techniques in Maxima.  Here’s a first approach:  conformal maps.  The code drawgridC below is limited in its versatility, but couldn’t be simpler to call.  Here are some examples,

/* plot the traces of complex valued function _f
 along grid lines in the square
 -1 <= Real(z) <= 1, -1<=Ima(z)<=1
 blue - images of _f along horizontal lines
 red - images of _f along vertical lines */
drawgridC(_f):=block(
[fx,fy,fxt,fyt,pl,j,ngrid],
fx:realpart(rectform(subst(z=x+%i*y,_f))),
fy:imagpart(rectform(subst(z=x+%i*y,_f))),
pl:[],
ngrid:19,
 for j:0 thru ngrid do(
 fxt:subst([x=-1+j*2/ngrid,y=t],fx),
 fyt:subst([x=-1+j*2/ngrid,y=t],fy),
 pl:cons([parametric(fxt,fyt,t,-1,1)] , pl)
 ),
pl:cons([color=red],pl),
 for j:0 thru ngrid do(
 fxt:subst([y=-1+j*2/ngrid,x=t],fx),
 fyt:subst([y=-1+j*2/ngrid,x=t],fy),
 pl:cons([parametric(fxt,fyt,t,-1,1)] , pl)
 ),
pl:cons([color=blue],pl),
pl:cons([xlabel="",ylabel=""],pl),
pl:cons([nticks=200],pl),
apply(wxdraw2d,pl)
);

Complex Factors of Real and Complex Polynomials – part 2

Last week I posted a Maxima utility to factor polynomials using complex roots.  The idea was to use the built-in solver solve to find each root a with multiplicity m, and then construct a factor of the form (z-a)^m.  Typically, solve can fail to identify some roots that are returned by the more robust to_poly_solve.  My problem at that time was that I didn’t know how to access the multiplicities for the roots returned by to_poly_solve.

I’ve written a utility to determine multiplicities for the roots returned by to_poly_solve. With that, here’s an updated version of factorC that takes advantage of the more robust solver.

/* factor a polynomial using all its complex roots */
factorC(_f,_z):=block(
[s,n,m,fp,j],
fp:1,
ss:to_poly_solve(_f,_z),
s:create_list(args(ss)[k][1],k,1,length(ss)),
m:tpsmult1(_f,_z,ss),
/*These lines were needed before I figured out how
 how to handle multiplicities with to_poly_solve
s:solve(_f,_z),
m:multiplicities,*/
n:length(s),
for j:1 thru n do 
 if lhs(s[j])#0 
 then fp:fp*(_z-(rhs(s[j])))^m[j],
 fp:fp*divide(_f,fp)[1],
fp
);

Roots, Multiplicities, and to_poly_solve

A few days ago I posted some code to factor and perform partial fraction decomposition using all available complex roots.  The approach is to find roots a — taking note of the multiplicity m of each —  and then construct factors of the form (z-a)^m.

As I’ve noted before, the built-in command solve can fail to find solutions.  The more robust Maxima solver is to_poly_solve.    That said, after a call to to_poly_solve, I don’t know how to get my hands on the multiplicities.  So, I’ve written a little utility to compute and return a list of multiplicities for the roots returned by to_poly_solve.  Here’s an example, followed by the code.  Notice that because to_poly_solve returns solutions using the %union construction, extracting the jth root requires something like root:rhs(part(s,j,1))

tpsmult1

/* tpsmult1 assumes results of to_poly_solve applied to _f, 
 are passed in the argument s */
tpsmult1(_f,_z,s):=block(
[d,root,nroots,mult,m1],
nroots:length(s),
mult:[],
for j:1 thru nroots do (
root:rhs(part(s,j,1)),
m1:rootmult(_f,_z,root),
mult: endcons(m1,mult)
),
mult
);

/* find the multiplicity of a root _a for the polynomial _f */
rootmult(_f,_z,_a):=block(
[d,val,m,mm],
m:1,
d:divide(_f,(_z-_a)),
val:subst(_z=_a,part(d,1)),
if (val=0) then (
 mm:rootmult(part(d,1),_z,_a),
 m:m+mm
 ),
m
);