Plotting with Maxima

Maxima is about computation, mainly symbolic computation. However, when you talk about math software, the first thing that comes to the mind of most people are beautifully rendered 3-D surfaces. The prototypical image is the sombrero-like surface you see below.


I've added it here just as an eye-catcher, because that's not where I want to begin. I'd like to start exploring the plot commands in Maxima with something simpler, but not as simple as the usual sines, cosines and polynomials. I wanted something a bit more challenging and at the same time related to a real problem. After searching for a while I came across Plank's Law, which every physics student studies in his or her first class of Quantum Mechanics.

Plank's Law gives the black-body radiation as a function of temperature and either frequency or wavelength. It can be expressed in a few different ways, depending on whether you want the density per solid angle or the total energy. This last form is commonly expressed as

You can find a good plot of this function for various temperatures in the Wikipedia article about Plank's Law, here.

So how do we go about plotting it with Maxima? There are basically two aspects to it. First you must encode the equation in a way that Maxima will understand. In other words, you must code it using Maxima's syntax. Second, you must use the plot2d command giving it proper arguments so that it can perform its job.

Instead of giving you the final result, I think it's more interesting if I tell you about some of my attempts and their gradual improvements. So we will do it in three steps.

Maxima functions and definitions can be kept in text files and loaded as needed with the "batch()" function. I keep a small batch file with the definitions of useful physical constants. In this example we will need

c: 2.9979e8$ /* Speed of light */
h: 6.6256e-34$ /* Plank's constant */
k: 1.3805e-23$ /* Boltzmann's constant */

You use ":" to separate a variable's name from its initial value. In this case the variables will actually be constants.

In a separate file I entered my first attempt:

/* First attempt at Planck's law */
/* Helper constants */
k1: 8 * %pi * h * c$
k2: (h * c) / k$

/* Function definition */
E(l,T) := k1 * (l ^ -5) * (1 / (exp(k2/(l * T)) - 1))$

/* Plot command */
plot2d(E(l,5500),[l,100e-9,2000e-9])$

Note that I defined two auxiliary constants, k1 and k2. Not only do they make the equation look simpler but they also avoid some unnecessary computation. Remember that Maxima will use the formula we give it to calculate the function's value in a range of very close points. So using these pre-calculated constants it won't need to calculate them over and over again for each point in the x-axis.

A few observations about the function definition. First note that it uses ":=" instead of ":" as we did for variables. Second, exponentiation is denoted by "^". It would have been nice to use the character lambda instead of the lowercase "l" and have wxMaxima display it as a Greek character, but I don't know how to do this. Maxima and wxMaxima do this nicely for other Greek letters, like alpha and theta, but somehow it does not work for lambda. Finally, operator precedence can be tricky so it's safer to use parentheses when in doubt.

As for the plot2d command, in its simpler form it takes two arguments. The first one is the expression we want to plot and the second one is the range of the independent variable. So we are asking it to plot the Energy at temperature T=5500K, for lambda varying from 100 to 2000 nm.

This was the resulting plot:


Not bad for a first try! Note that Maxima automatically scales the vertical axis so that the whole graph fits in, but the units for both axes are somewhat cryptic. In our second attempt we will try to display more friendly axes. Here's the improved code:

/* Second attempt at Planck's law */
/* Helper constants */
k1: 8 * %pi * h * c * 1e-3$
k2: (h * c) / k$

/* Function definition */
E(l,T) := k1 * ((l * 1e-9) ^ -5) *
(1 / (exp(k2/((l * 1e-9) * T)) - 1))$

/* Plot command */
plot2d(E(l,5500),
[l,100,2000],
[xlabel,"lambda [nm]"],
[ylabel,"E(lambda) [kJ/nm]"])$

In order to make the vertical axis look like the one in the graph from the Wikipedia page we have scaled the k1 constant by 1000. And since we want the horizontal axis to be in units of nm, we added the factor 1 x 10^-9 to each occurrence of lambda in the formula. Note the small quirk in Maxima's syntax. It uses the letter "e" to denote constants multiplied by a power of ten but uses a different symbol, "^", to denote exponentiation in expressions. No big deal, this is actually used in most (all?) programming languages, but it does look strange from a mathematical notation point of view.

Finally, we added some more arguments to the plot2d() command. They are the plotting options and take the form of lists enclosed by square brackets in which the first element is the option name. For this example we just attached labels to each of the axes, displaying their meaning and units. Here's the result:

We might say that we have accomplished our initial goal, at least for a single temperature. What remains to be done is add more temperatures, each in a different color. This will be very easy, since, for our convenience, the plot2d() commands accepts a list of expressions to plot as its first argument! So far we have given it a single expression, viz. E(l,5500), but now we will give it a list: [E(L,3500),E(L,4000),E(L,4500),E(L,5000),E(L,5500)]. Here I started using uppercase "L" for lambda because the lowercase version could be easily mistaken by the digit 1.

/* Final version of Planck's law plot */
/* Helper constants */
k1: 8 * %pi * h * c * 1e-3$
k2: (h * c) / k$

/* Function definition */
E(L,T) := k1 * ((L * 1e-9) ^ -5) *
(1 / (exp(k2/((L * 1e-9) * T)) - 1))$

/* Plot command */
plot2d([E(L,3500),E(L,4000),E(L,4500),E(L,5000),E(L,5500)],
[L,100,2000],
[xlabel,"lambda [nm]"],
[ylabel,"E(lambda) [kJ/nm]"],
[color,black,cyan,blue,green,red],
[legend,"T=3500K","T=4000K","T=4500K","T=5000K","T=5500K"]);

Besides, I have added more options to plot2d(). We now specify which colors are to be used for each plotted expression and which legend to attach to it. Unfortunately the legend does not appear near the plot but on the right-hand upper corner. So here's the final result:


Ok, so now I'm satisfied. This was a good exercise and the result looks very much like the plot from Wikipedia. But wait! E(L,T) is obviously a two-variable function although we have been using it like a single variable function by choosing only a few fixed values for the temperature. This is like a contour map of a 3D surface. We can do better. We can use Maxima's plot3d() command to have a continuous view of the Energy function for a range of temperatures! Here's how we do it:

plot3d(E(L,T),
[L,100,2000],
[T,1500,5500],
[zlabel,"E(L,T)"],
[legend,"Spectral Energy Density"],
[grid,50,50]);

It's really that simple! We just added a new axis label and the grid option to make the surface look smoother (although that was not necessary). Here's how it looks:


I had never seen Planck's function plotted like this. I guess the plotting software available around 1900, when Planck created it, was only capable of doing 2-D plots :-) This 3-D plot gave me a new insight into the problem. I immediately saw that the energy is not bounded by the temperature. When you look at the 2-D contour plot, perhaps you could overlook the fact that temperature can increase more and more and that the energy goes the same way. Not so in the 3-D plot, where temperature is one of the axes!

By the way, the sombrero plot of the beginning of this post was done with

plot3d(cos(x^2+y^2)*exp(-(x^2+y^2)/5),
[x,-4,4],
[y,-4,4],
[grid,100,100]);