Maxima Codes Fortran90 for Me

Maxima is great. In the magnetron post we integrated the differential equation describing an electron's motion in a couple of different ways. The implicit methods required a matrix inversion and then matrix-vector multiply for each time-step. Instead of writing out and coding the expressions describing that operation by hand (painful bugs quite likely, ouch!) it is much easier to let Maxima code the Fortran for us. If you have to ask "why Fortran?", then you are obviously some kind of quiche eating software weenie. If instead, you thought, "of course Fortran, what else would Maxima code?", please read on, you number crunching stud (or studette).

Start up a Maxima session (I run Fedora 9, but I'm sure most other distributions have it, I think they even have a windows installer). It's easy to define our matrix operator, invert it, and apply it:

L : matrix(
[1,0,0, -dt, 0, 0],
[0,1,0, 0, -dt, 0],
[0,0,1, 0, 0, -dt],
[0,0,0, 1, a*B[3]/c,-a*B[2]/c],
[0,0,0,-a*B[3]/c, 1, a*B[1]/c],
[0,0,0, a*B[2]/c,-a*B[1]/c, 1] );

Linv : invert(L);

X: [x,y,z,vx,vy,vz];
P : a*[0,0,0,E[1],E[2],E[3]];
RHS : X - P;

LinvRHS : factor( expand( Linv.RHS ) );


Simple as that we have a symbolic representation of our inverted operator acting on our state vector. Now we want Fortran90 expressions.

load("f90");

fname : "back_euler_inv_op.f90";
f : openw(fname);

printf( f, "! These fortran expressions generated automatically from
! the maxima batch file magnetron.mac ~%");

for i : 1 thru length(LinvRHS) do
( printf( f, "x_n1(~d)=",i),
with_stdout( f, f90(LinvRHS[i,1]) ) );

close(f);


There is a small bug in the "f90" module (as of version 5.15) that causes it to use the Fortran77 style continuation when writing out a matrix, so that is why each element must be written out with a separate call to f90().

Well, that's neat, but I need to appease the quiche eaters, so wouldn't it be better if these expressions, and our inner loops (the time integration loop for this case) in Fortran were callable from something flexible and trendy like Python? That's where F2Py comes in. Generating a module for Python is one simple call:

f2py -c back_euler_sub.f90 -m back_euler

Where my subroutine using the Maxima generated expressions is in the file "back_euler_sub.f90" and the module will be named "back_euler".

The python script to use our new module is quite short:

from numpy import *
from back_euler import *

B = zeros( 3 )
E = zeros( 3 )
x_n = zeros( 6 )
x_n1 = zeros( 6 )
q = -4.8032e-10
m = 9.1094e-28
c = 2.9979e10
dt = 3e-10
nt = 50
E[1] = -1.0/3.0
B[2] = 33.6412

X = zeros( (6,nt), dtype='double', order='F' )

X = back_euler(nt,dt,q,m,c,E,B)

Now our low-level number crunching and inner loop is in nice, fast compiled F90 and we can call it from a high-level, object-oriented language to do things like zero finding (as in the magnetron post), or optimization or something else entirely.