Month: February 2017

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**

linearalgebrapackage

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

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

Solving the matrix vector equation Ax=b in Maxima

*Upadate:  I’ve implemented a Maxima matrix-vector equation solver with simpler Maxima-specific algorithm in a later post.  That method is based on the built-in function linsolve().  In that post I show an example that breaks linsolve() but that is handled correctly by the backsolve() method.

Is there really not a solver in Maxima that takes matrix A and vector b and returns the solution of Ax=b ?  Of course we could do invert(A).b, but that ignores consistent systems where A isn’t invertible…or even isn’t square.

Here’s a little function matsolve(A,b)  that solves Ax=b for general A using the built-in Gaussian Elimination routine echelon(), with the addition of a homemade backsolve() function.  The function in turn relies on a little pivot column detector pivot() and my matrix dimension utility matsize(). This should include the possibilities of non-square A, non-invertible A, and treats the case of non-unique solutions in a more or less systematic way.

matsolve

matsolve(A,b):=block(
   [AugU],
   AugU:echelon(addcol(A,b)),
   backsolve(AugU)
 );

backsolve(augU):=block(
 [i,j,m,n,b,x,klist,k,np,nosoln:false],
 [m,n]:matsize(augU),
 b:col(augU,n),
 klist:makelist(concat('%k,i),i,1,n-1),
 k:0,
 x:transpose(matrix(klist)),
 for i:m thru 1 step -1 do (
   np:pivot(row(augU,i)), 
   if is(equal(np,n)) then
     (nosoln:true,return())
   else if not(is(equal(np,0))) then
     (x[np]:b[i],
     for j:np+1 thru n-1 do
       x[np]:x[np]-augU[i,j]*x[j])
    ),
 if nosoln then 
    return([])
 else 
   return(expand(x)) 
)$

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

pivot(rr):=block([i,rlen],
 p:0,
 rlen:length(transpose(rr)),
 for i:1 thru rlen do( 
 if is(equal(part(rr,1,i),1)) then (p:i,return())), 
 return(p)
)$

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

eigs

For a matrix with a full set of eigenvectors but eigenvalues of multiplicity greater than one, the lines above fail.  A version of the above that works correctly in that case could look like:

eigdiag(v):=apply('diag_matrix,flatten(makelist(makelist(part(v,1,1,j),i,1,part(v,1,2,j)),j,1,length(part(v,1,1)))));

eigU(v):=transpose(apply('matrix,makelist(makelist(flatten(part(v,2))[i],i,lsum(i,i,part(v,1,2))*j+1,lsum(i,i,part(v,1,2))*j+lsum(i,i,part(v,1,2))),j,0,lsum(i,i,part(v,1,2))-1)));

BDF2a: An adaptive-stepsize Stiff ODE solver in Maxima

**March 2021 update** Maxima has a full-featured stiff solver called Eulix now, thanks to Helmut Jarausch.

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

Since 2011, Maxima has included the user-contributed numerical ODE solver rkf45() created by Panagiotis Papasotiriou.  This implementation of the fourth (and fifth) order Runge-Kutta-Fehlberg embedded method features adaptive timestep selection and a nicely optimized function evaluation to make it run pretty fast in Maxima.  As an explicit RK method, it is suitable for non-stiff equations.  To use it in Maxima, you must first load(rkf45).

For years in my differential equations class, we’ve test-driven the various numerical ODE methods in MATLAB and explored their behavior on stiff systems.  This year we’ve switched from MATLAB to Maxima, and I’ve written a stiff solver for that purpose.  It’s an adaptive stepsize  implementation of the second-order Backward Difference Formula.  The BDF.mac file can be found here.  The package contains a fixed-stepsize implementation BDF2()  and the adaptive stepsize BDF2a().

I’ve tried to make the calling sequence and options as close to rkf45() as possible.  I have not made any attempt to optimize the performance of the method.  BDF methods are implicit, meaning that a nonlinear equation (or system of equations) must be solved at each step.  For that I use the user-contributed Maxima function mnewton().

Additionally, note that I included my homemade help utility so that once BDF.mac is loaded, help(BDF2) and help(BDF2a) display help lines in the console window. For simplicity in comparing the methods, the load(rkf45) command is included in BDF.mac, as is an rkf45 help file for use with help().   Finally, I included in the package my little solution plotting utility wxtimeplot() which takes the “list of lists” output of BDF2(), BDF2a() and rkf45() and plots the dependent variables vs the independent variable along with a log-scale plot of the stepsize.

Below is an example that illustrates how the the BDF method performs relative to rkf45 for an increasingly stiff problem.

The equation is

 \dot{y} = -\lambda (y-\sin (10t)-t)+10 \cos (10t)+1 .    

For small values of  \lambda , the explicit rkf method’s adaptive stepsizes result in much more efficient solutions for a given accuracy level.  However for large \lambda the implicit BDF method is comparatively much more efficient, in terms of number of steps, and  number of function evaluations, although in its current implementation, the BDF method takes more time.

Here’s the results for \lambda=50:

stiff50
stiff50plot

But for \lambda=5000, the performance situation is reversed, with the BDF formula essentially immune from the effects of extra stiffness, while the rkf method requires an order of magnitude more steps.

stiff5000
stiff5000plot

References:

E. Alberdi Celaya , J. J. Anza Aguirrezabala , and P. Chatzipantelidis.  “Implementation of an adaptive BDF2 formula and comparison with the MATLAB Ode15s”,  Procedia Computer Science,  December 2014

Using Maxima To Do Some Math Described in the Movie “Hidden Figures”

I looked in the NASA archives for the work of Katherine Johnson described so dramatically in the film “Hidden Figures”.  I found this 1960 paper:

nasa_title

In particular, I was struck by the film’s use of “Euler’s Method” as the critical plot point, and sadly couldn’t find anything about the use of that method in the archive, nor mention of it in the text of  Margot Lee Shetterly’s book on which the film “Hidden Figures” is based.

I did however read in the paper about an “iterative procedure” needed to solve for one variable.

Turns out the sequence of equations described in the report, when assembled into an expression in the single variable of interest, defines  a contraction mapping and so converges by functional iteration, as is numerically demonstrated in the paper.

At the bottom I’ve linked a fuller treatment in a .wxmx file.  Briefly, the authors describe the problem:

nasaeq19

nasaeq20

nasaeq89

nasa_iteration

Here’s everything from the paper needed to do that in Maxima:

nasaconstants1

nasaconstants2

nasaconstants3

I put together equations 19,20,8, and 9, together with all those constants defined by the paper and customized trig function that act on angles measured in degrees, into an expression of the form

  \theta_{2e}=F(\theta_{2e})

and applied the built-in Maxima numerical solver find_root().  The result agrees reasonably well with the result  \theta_{2e}=50.806 quoted in the paper.

nasaiteration_findroots

In the above, I computed the derivative of the mapping F evaluated at the solution and found that the value is much less than one, showing it is a contraction.

Here’s my full wxMaxima session.

Maxima wins Sourceforge “Community Choice Project of the Month”

 

With over 300,000 downloads per year, Maxima is gaining popularity as a Computer Algebra System comparable to commercial systems like Mathematica and Maple.

The Maxima project was recognized with a community choice award at Sourceforge for the month of February.   Previously, Maxima was chosen as “Staff Pick Project of the Month” in November, 2015

https://sourceforge.net/blog/february-2017-community-choice-project-of-the-month-maxima/

trigexpand(), trigsimp() and an oscillatory solution of a 2nd order linear differential equation

I was working on a differential equations homework problem and it took me a few tries to remember how to simplify a trig expression.  So, I thought I’d put the results here to remind myself and others about trigexpand().

The second order linear equation is

 y''+y'+y=0 \;\;\;\;\;\; y(0)=1, y'(0)=0

odesol

Now I asked my students to use the identities

  A= \sqrt{C_1^2 + C_2^2}  and  \tan \phi=C_2/C_1

to write the equation equivalently as

 y=A e^{-x/2} \cos \left( \frac{\sqrt{3}}{2}  x -\phi \right)

odeoscillparams

Finally, I wanted to use Maxima to show these two expressions are equivalent.  I reflexively used trigsimp() on the difference between the two expression   (the name says it all, am I right?)  but the result was complicated in a way that I couldn’t recover…certainly not zero as I hoped.  A few minutes of thought and I remembered trigexpand() which, followed by a call to trigsimp(), gives the proof:

odetrigexpand