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

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
-



## Jordan Canonical Form in Maxima

After not easily finding such a thing from a cursory search of the Maxima documentation, I spent a few hours over the weekend beginning to write a Maxima function to compute, for any given square matrix $M$, an invertible matrix $P$ so that

$P^{-1}MP = J$

where $J$ is the Jordan matrix that displays the eigenvalue/vector structure of $M$.

It took several searches for me to find, but of course there’s already such a  function — with a not so easily searched-for name — in the diag package:   ModeMatrix()

To see just the matrix $J$, diag provides jordan() and dispJordan()

## Solve Ax=b in Maxima, part 2

In a previous post, I included my little coding project to implement a general backsolve() function to use with the built-in maxima matrix function echelon(), producing an easy-to-call matrix solver matsolve(A,b).  The result is meant to solve a general matrix vector equation $Ax=b$ , including cases when $A$ is non-square and/or non-invertible.

Here’s a quicker approach — convert the matrix into an explicit system of equations using a vector of dummy variables, feed the result into the built-in Maxima function linsolve(), and then extract the right hand sides of the resulting solutions and put them into a column vector.

The two methods often behave identically, but here’s an example that breaks the linsolve() method, where the backsolve() method gives a correct solution:

*Note, I’ve found that the symbol rhs is a very popular thing for users to call their problem-specific vectors or functions.  Maxima’s  “all symbols are global” bug/feature generally wouldn’t cause a problem with a function call to rhs(), but the function map(rhs, list of equations)  ignores that rhs() is a function and uses user-defined rhs.  For that reason I protect that name in the block declarations so that rhs() works as expected in the map() line at the bottom.  I think I could have done the same thing with a quote: map(‘rhs, list of equations).

matsolve2(A,b):=block(
[rhs,inp,sol,Ax,m,n,vars],
[m,n]:[length(A),length(transpose(A))],
vars:makelist(xx[i],i,1,n,1),
Ax:A.vars,
inp:makelist(part(Ax,i,1)=b[i],i,1,n,1),
sol:linsolve(inp,vars),
expand(transpose(matrix(map(rhs,sol))))
);

## A Little Maxima Function to Find the Dimensions of a Matrix

**Update**  I didn’t find it it in documentation for quite a while, but there is a built-in Maxima function matrix_size()  in the package linearalgebra that does what this little one-liner does**

I really wanted a Maxima function that works something like MATLAB size() to easily determine the number of rows and columns for a matrix $M.$  In Maxima, length(M) gives the number of rows, and so length(transpose(M)) gives the number of columns.  I put those together in a little widget matsize() that returns the list [m,n] for an $m \times n$ matrix $M:$

matsize(A):=[length(A),length(transpose(A))];

## Extracting Matrices from the Output of eigenvectors()

In class I sometimes need to use matrices of eigenvalues and eigenvectors, but the output of eigenvectors() isn’t particular helpful for that out of the box.

Here are two one-liners that work in the case of simple eigenvalues.  I’ll post updates as needed:

First eigU(), takes the output of eigenvectors() and returns matrix of eigenvectors:

eigU(v):=transpose(apply(‘matrix,makelist(part(v,2,i,1),i,1,length(part(v,2)),1)));

And eigdiag(), which takes the output of eigenvectors() and returns diagonal matrix of eigenvalues:

eigdiag(v):=apply(‘diag_matrix,part(v,1,1));