Using Maxima Output in Octave

Introduction

This is the third in a series of articles on Octave starting with Octave, An Alternative to the High Cost of MATLAB. Octave is a free, both free as in beer and free as in speech, MATLAB compatible numerical programming tool available under the GNU General Public License. In part because MATLAB has become the de facto industry standard for numerical programming, Octave is of particular interest to individuals, companies, and organizations engaged in numerical and mathematical programming and research and development.

Octave has some limitations. The base Octave tool has no symbolic manipulation features. It is not a computer algebra system (CAS) such as Mathematica or Maple. It cannot, for example, perform symbolic integration, symbolic differentiation, factor polynomials, and so forth. Octave does have the symbolic toolbox available through the Octave Forge repository of Octave toolboxes, but the symbolic toolbox is quite limited. A better option for symbolic manipulation tasks is to use the Maxima computer algebra system in combination with Octave. Octave also lacks the ability to generate TeX or LaTex mathematical output, the de facto standard of mathematical publication. Maxima can also generate TeX output for inclusion in papers or WordPress blog posts.

Maxima

Maxima is a computer algebra system descended from MACSYMA, one of the original computer algebra systems. MACSYMA was developed at MIT in part for use in theoretical physics. Maxima is available both as source code and pre-compiled binaries for all three major computer platforms: Unix/Linux, Microsoft Windows, and Mac OS. Maxima is free software, both free as in beer and free as in speech, available under the GNU General Public License (GPL). wxMaxima is a Graphical User Interface (GUI) for Maxima available as both source code and pre-compiled binaries for all three major computer platforms. wxMaxima has human readable menu items and buttons for many common symbolic manipulation and mathematical functions. wxMaxima also has “notebooks” similar to Mathematica notebooks. There is considerable documentation on Maxima; interested readers are referred to the excellent online and published documentation on Maxima. This article is focused on using Maxima as an adjunct to Octave.

wxMaxima on a Macintosh

Maxima can perform both symbolic differentiation and integration. Symbolic differentiation is illustrated in the screen shot of Maxima above. Some optimization algorithms, used, for example, for model fitting, require the derivative of the function being optimized; the function is usually being minimized. If the function is rather complex, deriving the derivatives of the function with respect to the parameters over which the optimization is performed by hand can be time consuming, tedious, and error prone. The author has used Maxima successfully to perform the differentiation of a model function. One can then convert the Maxima output, the derivative produced by Maxima’s symbolic differentiation, into an expression that can be used in Octave by using the fortran(expression) command in Maxima. The Maxima fortran command generates FORTRAN code for the Maxima expression. In many cases, the FORTRAN expressions are identical to Octave/MATLAB mathematical expressions. In some cases, the FORTRAN code generated by Maxima must be edited slightly to create valid Octave/MATLAB code.

For example, the Cauchy-Lorentz distribution is a commonly used mathematical model of a peak.

The Cauchy-Lorentz distribution is the frequency response of a forced-damped harmonic oscillator. It is widely used in physics, mathematics, and engineering under a number of different names.

Cauchy Lorentz Distribution

In using the Cauchy Lorentz distribution to model a peak in data, one typically wants to determine the values of the parameters A, mu, and W representing the magnitude of the peak (A), the position of the peak (mu), and the width of the peak (W) that best fits the data. To do this, some model fitting algorithms need the derivatives of the Cauchy Lorentz distribution with respect to each parameter A, mu, and W.

Differentiation of Cauchy Lorentz in Maxima

This is the derivative of the Cauchy Lorentz distribution with respect to the width parameter W from symbolic differentiation in Maxima:

This derivative is moderately complex. Calculating this derivative by hand is time consuming and error prone. Imagine computing the derivative of an extremely complex mathematical model with hundreds of terms by hand. The probability of error even by a highly-skilled mathematician is very high. It was for this reason that tools like Maxima were developed.

This is the FORTAN code generated by applying the Maxima fortran(expression) function to the Maxima expression for the derivative of the Cauchy Lorentz with respect to the width W

2*(x-mu)**2*A/(((x-mu)**2/W**2+1)**2*W**3)

This is actually valid Octave/MATLAB code. If the variables x, mu, A, and W are scalar variables in Octave, this FORTRAN expression will evaluate correctly. Here is the calculation in Octave when x, mu, A, and W are all scalar variables with the value 1.0.

octave-3.2.4.exe:25> x = 1
x = 1
octave-3.2.4.exe:26> 2*(x-mu)**2*A/(((x-mu)**2/W**2+1)**2*W**3)
ans = 0
octave-3.2.4.exe:27>

However, in Octave the variable x is often a vector. If x is a vector, the expression above will produce an error in Octave:

octave-3.2.4.exe:25> 2*(x-mu)**2*A/(((x-mu)**2/W**2+1)**2*W**3)
error: for A^b, A must be square
octave-3.2.4.exe:25>

The reason for this error is that in Octave and MATLAB some of the operators such as * and / are not by default interpreted as element by element operators when applied to vectors and matrices. For example, the operator * is matrix multiplication by default in Octave and MATLAB. An element by element operator is an operator that is applied separately to each element in each vector or matrix that is an operand. In Octave and MATLAB, the element by element operators are .*, ./, .+, .-, and so forth. For example, if one has two vectors a and b in Octave, the operator * will give an error:

octave-3.2.4.exe:34> a
a =

1  2  3  4  5

octave-3.2.4.exe:35> b
b =

1  2  3  4  5

octave-3.2.4.exe:36> a * b
error: operator *: nonconformant arguments (op1 is 1x5, op2 is 1x5)
octave-3.2.4.exe:36>

However, in Octave and MATLAB, one can multiply each element of each vector by the corresponding element of the other vector using the element by element operator .* thus:

octave-3.2.4.exe:36> a .* b
ans =

1   4   9  16  25

octave-3.2.4.exe:37>

In the output of the element by element (or elementwise) operator .*, the first element is 1*1, the second element is 2*2, and so forth.

Thus, the FORTRAN expressions generated by Maxima are not valid Octave/MATLAB code for vectors and matrices, only for scalar variables. One can convert the FORTRAN expression to a valid Octave expression for vectors by converting the non-elementwise operators to element by element operators where the operands are vectors, usually by adding a preceding dot. Here is the edited code for the example derivative:

function [result] = mydiff_edited(x, A, mu, W)
% function [result] = mydiff(x, A, mu, W)
% FORTRAN code from wxMaxima edited to support an array x as input
%

result = 2*(x-mu).**2*A ./ (((x-mu).**2 ./ W.**2+1).**2*W**3)

end

If x is defined as a vector in Octave such as:

x = [0.0:0.1:10.0]; % x values from 0.0 to 10.0 in steps of 0.1

and then this function is used to compute the value of the derivative of the Cauchy Lorentz distribution with respect to the width parameter W:

data = mydiff_edited(x, 1.0, 1.0, 1.0);
plot(data);

Plot of Vector Derivative

Maxima can also generate valid TeX code for mathematical publication through its built in tex(expression) command. This command can also be invoked through a menu item in the wxMaxima GUI (shown below).

wxMaxima with TeX Menu Item Displayed

Some TeX generated by Maxima:

tex(1/(1+x^2));

generates the TeX code:

$${{1}\over{x^2+1}}$$

which displays in WordPress after removing the $$ tags which WordPress does not need as:

All of the mathematical formulas in this article are TeX generated by Maxima in this way.

Conclusion

Octave has extensive numerical analysis and programming features. Octave has the special advantage that it is mostly compatible with MATLAB, which is currently the de facto industry standard for numerical programming. Most Octave scripts will run under MATLAB and many MATLAB scripts will run under Octave with no changes. If a user or software developer has an occasional need for symbolic manipulation features such as symbolic integration and differentiation, one can use Maxima as an adjunct to Octave. Similarly, one can use Maxima to generate TeX code for mathematical publications. If one needs to perform extensive symbolic manipulation, one may need to use Maxima or similar tools as one’s primary tool.

Both Octave and Maxima have the advantage that they are free, both free as in beer and as in speech, and available as source code. There are many cases where a merger, change in corporate strategy, bankruptcy, or even the whim of an executive has resulted in a proprietary development platform being discarded or deprecated to the detriment of end users, developers, and other customers. A well known example is FoxPro, once one of the leading database programs, which Microsoft acquired and has now announced will be discontinued in favor of Microsoft’s other database products. In contrast, open source development tools such as Octave and Maxima can be kept alive and indeed improved by their end users, developers, and customers if needed.

© 2011 John F. McGowan

About the Author

John F. McGowan, Ph.D. is a software developer, research scientist, and consultant. He works primarily in the area of complex algorithms that embody advanced mathematical and logical concepts, including speech recognition and video compression technologies. He has extensive experience developing software in C, C++, Visual Basic, Mathematica, MATLAB, and many other programming languages. He is probably best known for his AVI Overview, an Internet FAQ (Frequently Asked Questions) on the Microsoft AVI (Audio Video Interleave) file format. He has worked as a contractor at NASA Ames Research Center involved in the research and development of image and video processing algorithms and technology. He has published articles on the origin and evolution of life, the exploration of Mars (anticipating the discovery of methane on Mars), and cheap access to space. He has a Ph.D. in physics from the University of Illinois at Urbana-Champaign and a B.S. in physics from the California Institute of Technology (Caltech). He can be reached at [email protected].

Sponsor’s message: Check out Math Better Explained, an insightful ebook and screencast series that will help you see math in a new light and experience more of those awesome “aha!” moments when ideas suddenly click.