## Partial Fraction Decomposition Using Complex Factors

A few days ago I wrote about my solution for factoring a polynomial into factors corresponding to complex roots,  factorC.

One of the uses for that utility is to perform a partial fractions decomposition of rational functions.  Here’s a complex partial fraction function that uses factorC and the built-in Maxima function partfrac.

partfracC(_f,_z):=block(
[d,fd],
d:denom(_f),
fd:factorC(d,_z),
partfrac(1/fd,_z)
);

## Complex Factors of Real and Complex Polynomials

I’m getting ready for my fall Complex Variables class.  I noticed that the built-in Maxima function residue doesn’t reliably do the right thing.  My goal is to make some improvements to  Maxima residue calculations in Maxima over the course of the next month.

As I started to look at some test cases, I realized I didn’t know how to factor a polynomial into complex factors.  In the simplest case:

but I wanted to see

Maybe someone will find this and let me know that there’s a simple way to make that happen using existing Maxima commands.  Until then, I’ve written a little utility  to identify the roots (both real and complex) of a polynomial and return a factorization.  Here’s an example, and the code.  Notice that it only works as well as the root finder solve.  I tried to upgrade to the more robust to_poly_solve,  but I don’t yet know how to handle multiplicities in that case.

factorC(_f,_z):=block(
[s,n,m,fp,j],
fp:1,
/* This commented code was meant to use the
more robust solver to_poly_solve, but
I couldn't understand how to handle multiplicities
ss:args(to_poly_solve(_f,_z)),
s:create_list(ss[k][1],k,1,length(ss)),*/
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
);

## Computing and Visualizing Path Integrals of Vector Fields

A few months ago, I posted my path integral functions, which are included in the MATH214 package.  Recently, I came across something I’d been looking for:  a Maxima utility for visualizing vector fields.  Its in the Maxima/Share directory under drawutils.

Written in 2010 by Donald J Bindner, the commands plot_vector_field and plot_vector_field3d do almost everything I was looking for.  The drawback is that I wanted to plot the vector fields along with the integration path.  I modified those two commands slightly into versions called make_vector_field  and make_vector_field3d to produce the lists of vectors for plugging into draw2d and draw3d, so that I could include the vector fields in  bigger graphics calls.  My modifications are available here.  The package includes my home-baked help utility.

Here’s what the path integral command and the vector field generator look like on an example from chapter 13 of Stewart’s “Calculus Concepts and Contexts”

## Solvers!

A few months ago, when I was very new to Maxima, I posted my solutions to several multivariable optimization problems, and included a standard constrained optimization example for which the built in command solve failed.

Here’s the problem:

In my multivariable calculus class, we’ve just encountered that problem again.  What my students and I discovered was that the order of the variables as specified by the user in the argument list caused the command to either succeed or fail, as you can see from the first two lines below:

I wrote to the experts at maxima-discuss@lists.sourceforge.net.  They replied en masse very quickly and were very generous with their time and advice. They verified the problem, pointing out that solving systems of equations is really complicated, and that the solve program seeks to solve for the variables one at a time and things can go wrong if an inauspicious order of variables is followed.

The experts also tried our problem with the package to_poly_solve, and happily it succeeded without the hassle of needing to specify the variables in the right order.  For solve users who reference the solutions it produces as elements of a list, direct use of to_poly_solve presents a  challenge due to output in a format different from that of solve  — t0_poly_solve returns solutions as a set with %union.

I’ve written a little wrapper Ksolve.mac that is called like solve as illustrated above.

In Ksolve if an initial call to solve fails, the process automatically upgrades to to_poly_solve in several variations, and  (hopefully) returns solutions in the same format as solve.

## Maxima for maxima (and minima) of functions z=f(x,y)

Together with the 3D graphing capabilities of Maxima, we can bring  symbolic differentiation and the numerical solver to bear when we seek local extrema of a surface.

Here is a link to the html export of a wxMaxima session where I work on two examples from my multivariable calculus class.  And here is the wxMaxima session.

## The MATH214 package for multivariable calculus

I’ve written a suite of Maxima commands for use in multivariable calculus class.  The package includes:

cross, dot, len,
unitT, unitN, unitB, curvature,
integrate2, integrate3,
integratePaths, integratePathv2, integratePathv3,
integrateSurf,
grad, div, curl

To use these commands:

2. extract (unzip) the contents into a directory you can find later. (I think macOS does this step behind the scenes, but in windows you’ll need right-click and “Extract All…”)
3. In wxMaxima, select File—Load Package…  then navigate to the directory in step 2. above and select MATH214.mac   The result should be an automatically generated input line similar to:

Included in that package is my home-baked help utility.  After loading, help(MATH214) returns a list of functions in the package, and help(<function_name>) returns a description and usage example for any of the functions named in the list above.

Examples showing these commands being executed in wxMaxima can be found here.  In addition, the .zip file also contains these examples in the  wxMaxima session file math214_testfunctions.wxmx that you can load into your wxMaxima session.

It is only necessary to download and extract MATH214.zip once, but you will need to load the package (step 3 above) in each new wxMaxima session you’d like to use any of these commands.

Note that these commands duplicate existing Maxima functionality (like dot(x,y) and x.y) or perform similarly to other packages (like vect).  The purpose here is for the commands in MATH214 to have a calling syntax that closely follows the way we have defined these operations analytically in class, while avoiding the unfortunate namespace  conflicts between the existing packages vect and draw that has been recently documented.

## 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:

And here’s another famous phase curve—the solution of the Brusselator:

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!

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:

## Path Integrals

For my multivariable calculus class, I wanted an easy-to-call suite of symbolic integrators for path integrals of the form

$\int_C f(x,y) dr$,

$\int_C {\bf F}(x,y)\cdot d{\bf r} = \int_C \langle P(x,y), Q(x,y) \rangle\cdot \langle dx, dy \rangle$, or

$\int_C {\bf F}(x,y,z)\cdot d{\bf r} = \int_C \langle P(x,y,z), Q(x,y,z), R(x,y,z) \rangle\cdot \langle dx, dy, dz \rangle$.

My overarching design idea was that the input arguments needed to look the way they do when I teach the course:

• a scalar field $f(x,y):R^2 \rightarrow R$ or a vector field  ${\bf F}(x,y): R^2 \rightarrow R^2$ or ${\bf F}(x,y,z): R^3 \rightarrow R^3$
• a curve $C$ defined by a vector-valued function ${\bf r}(t): R \rightarrow R^n, a\le t \le b$ where $n=2,3$ as appropriate.

It took me a while to work out how to evaluate the integrand along the path within my function.  Things that worked fine on the command line failed when embedded into a batch file to which I passed functions as arguments.  I ended up using subst, one variable at a time.  I’d like to be able to do this in a single command which can detect whether we’re in 2 or 3 dimensions so that I don’t need separate commands.

For now, here’s what I came up with along with some illustrative examples taken from Paul’s online math notes, that show how to call these new commands I, I2 and I3.

/* path integral of a scalar integrand f(x,y) on path r(t) in R^2, t from a to b */
I(f,r,t,a,b):=block(
[f1,f2,dr,Iout],
f1:subst(x=r[1],f),
f2:subst(y=r[2],f1),
dr:sqrt(diff(r,t).diff(r,t)),
Iout: integrate(f2*dr,t,a,b),
Iout
);
/* path integral of a vector integrand F(x,y) on path r(t) in R^2, t from a to b */
I2(H,r,t,a,b):=block(
[H1,H2,I],
H1:subst(x=r[1],H),
H2:subst(y=r[2],H1),
I: integrate( H2.diff(r,t),t,a,b),
I
);
/* path integral of a vector integrand F(x,y,z) on path r(t) in R^3, t from a to b */
I3(H,r,t,a,b):=block(
[H1,H2,H3,I],
H1:subst(x=r[1],H),
H2:subst(y=r[2],H1),
H3:subst(z=r[3],H2),
I: integrate( H3.diff(r,t),t,a,b),
I
);

Here’s an update:  a related maxima function for evaluating a complex integral

$\int_\Gamma f(z) dz$

where $f: C \rightarrow C$ and the curve $\Gamma$ is given by $r: R \rightarrow C$.

/* path integral of a complex integrand f(z): C --> C, on path z(t): R --> C, t from a to b */
IC(f,r,t,a,b):=block(
[f1,dz,Iout],
f1:subst(z=r,f),
dz:diff(r,t),
Iout: integrate(f1*dz,t,a,b),
Iout
);

## :=, ”(), define and div, grad, curl

I recently posted about : and :=  for defining functional expressions.  I’m starting to enjoy these emoji-like constructions 😉

This is another  colon-equals post.  This time for defining functions involving the maxima differentiation command diff.

Notice below that if we define a function with :=, the naive use of :=diff doesn’t produce a derivative with the expected results upon evaluation.

In fact, it’s a good thing that :=diff works like that.  The error with fp(3) above comes from the fact that we’ve actually defined an operator that differentiates the function with respect to the argument we pass…in the case above, differentiating with respect to the symbol u makes sense, while differentiating with respect to the constant 3 doesn’t.

So how to make the derivative function do what we want?  Two ways, that are subtly different, in ways I’m not completely sure of.  More about that when I learn more :-).

First is define,

Also you can use    ”()       quote-quote with parens around the whole right hand side:

I used define to write functions for vector valued 3D curves in an earlier post.   In figuring this out, I also learned that the :=diff form is really useful.  Below are three little functions in which I use :=diff to define the vector calculus operators grad, div and curl.  Notice that we pass the function f as an argument, and the :=diff form allows Maxima to differentiate them behind the scenes and return the results of the grad, div, and curl operators as you’d expect. These versions of div, grad and curl behave differently, and for me more as expected, than the functions of those names included in the Maxima vect package.  You can download the .mac file here.

/* Three Maxima functions for the multivariable calculus operators  grad, div, and curl
TheMaximaList.org, 2016
*/

grad(f,x,y,z):=[diff(f,x),diff(f,y),diff(f,z)]$div(f,x,y,z):=diff(f[1],x)+diff(f[2],y)+diff(f[3],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) ]\$

Here is a screenshot showing how to call these functions: