A tour of Maxima (with wxMaxima)

José A. Vallejo
Faculty of Science
State University of San Luis Potosí (México)
email: jvallejo@fc.uaslp.mx
URL: http://galia.fc.uaslp.mx/~jvallejo
Date: May 19 2019

ABSTRACT: In this document we try to illustrate the possibilities offered by the CAS (Computer Algebra System) Maxima, along the lines of the chapters "A Tour of Mathematica" that appeared in the printed versions of The Mathematica Book by Wolfram Research up to its 5th edition (2003). Although it is free software (as in 'freedom' and as in 'free beer'), the quality of Maxima is comparable to that of commercial programs such as Mathematica (TM) or Maple (TM). Its development began at the MIT (Massachusetts Institute of Technology) in the 60s, and it is nowadays active thanks to the dedicated work of a wide group of people around the world. Maxima admits many graphical interfaces (originally, it is invoked from the shell). The one we will use in this document is wxMaxima, although almost all of the results are independent of the interface (the only exception being the animations). This document has been created on a Linux platform running Slackware 14.2, with Maxima version 5.38.0 compiled against SBCL 1.3.5, and wxMaxima version 16.04.2.

REMARK: If the mathematical symbols are displayed in a sub-optimal way, please change the MathJax renderer of your browser: just click with the right button of the mouse on any of the equations below and, from the menu, select Math Settings - Math Renderer - SVG (for instance). The displayed output depends on which sets of fonts do you have installed in your system.

 1 Basic Arithmetics

Arithmetics can be done in Maxima as in a calculator. For instance, we enter 13+11 and Maxima returns the result, 24.

However, in contrast to a calculator, Maxima can work with exact computations. Here we give the value of 3^100. The circumflex ^ is the symbol for the exponentiation.

The 'float' function gives approximate numerical results (by default, with 15 digits precision). The percent sign % is used to make reference to the last output (older results can be accesed trough %th(n)):

Numerical results can be evaluated with any desired precision. The next command computes sqrt(10) with 50 digits, using "bfloat" (bigfloat number):

Maxima can also work with complex numbers. Here we show the rectangular form of (3+4i)^10 computed in two alternative forms. The complex unit i is denoted %i (other important constants implemented by default are %e, %pi or %phi, notice that all of them begin with a % sign).

However, expansion is not done by default: in Maxima (3+4i)^10 is an expression in its own right!

More operations with complex numbers. The syntax is self-explaining, let us only remark that the colon ':' is used to assign values to variables (for instance, in the first line we assign the label z1 to the complex number 2+3i).

As we are working symbolically, the result is not automatically simplified, we must explicitly ask for the simplification, in this case, to pass to the rectangular form:

Of course, it is possible to go back and forth between the rectangular and polar forms:

In these commands we have assigned a value to a variable through ':' as in z1:2+3*%i. Once that assignment is made, each time the symbol z1 appears Maxima will understand that it has the assigned value:

If we want to remove the assignment, we can use 'remvalue' (another, more drastic, option is to use 'kill')

Lists are very basic constructs in Maxima. They are collections of elements enclosed between square brackets, and they can be assigned to a variable name like A, mylist, etc. For instance (the dollar sign $ is used here to prevent Maxima from writing the output of a command, which in this case is just a copy of the input):

Operations with lists are immediate:

Also, most of the built-in functions of Maxima can operate directly on lists:

In many cases, it is necessary that the length of the lists be the same. Otherwise, Maxima will complain, as occurs if we try [1,2,3,4,5]*listB.

Functions of the argument x can be defined with the same syntax used in the blackboard, using ':=' as in the following example:

and then can be evaluated:

But Maxima offers another possibility for defining functions, called lambda functions, which are extremely powerful. The idea is to work with unnamed functions, as when one says 'consider the function that sends the variable argument u to the expression u^2+1'. This is done as follows:

Obviously, this can't be used as it stands, but we can say something like 'map the function that squares its argument and sums 1 to the result' to every element in the list [1,2,3,4], like this:

Isn't it magic?

Maxima can compute the numerical values of all the standard functions. As an example, consider the value at 14.5 of the Bessel function of the first kind J0:

Next, we are going to compute a real root of J0(x) near x=14.5, with a 20-digits precision. We will use the 'mnewton' routine for this (as its name says, it uses the Newton method). As this routine is not loaded by default, we call it first. Also, notice the syntax: usually when we pass an argument to a function that works with equations, an expression 'expr' is taken as the equation 'expr=0'

Some functions, as 'mnewton', return their output as a list (in this case, a list with only one element which is in turn a list, [x=14.93091770848779]. That is called  nested list). If we want to extract the result, we can do the following (another possibility is to use %[1](%) below, that is, [i] selects the i-th element of the list):

The reason for giving the output in the form of a list is that the command 'mnewton' has lots of options, some of them returning more onformation than just a value. Here, we solve a nonlinear system of equations:

Let us see the result:

 2 Elementary Algebra and Calculus

Before we begin, let's clear all the variable assignments we have done up to now:

We can check that no assignment is stored by looking at the output of 'values', which keeps track of those of the form a:something

and by looking at the output of 'functions', which keeps track of the definitions of the form f(x):=

Recall that Maxima can work not only with numbers, but also with algebraic expressions. Here we consider the expression 9(2+x)(x+y)+(x+y)^2:

The next command takes this expression, represented by %, raises it to the power of 3 and expands the result:

By factoring the output we obtain a much simpler form:

Factoring can be done in different numerical domains. By default, it is carried on inside the integers ring Z:

The following modification of this command executes the factorization in the Gaussian integers Z[i]:

It is possible to factor within other domains, it is enough to adjoin the appropriate roots (this is done by declaring the minimal polynomials that they satisfy). For instance, to factor in Z[sqrt(5)], we would do

Due to their importance in many topics, Maxima has a good deal of functions related to polynomials:

The coefficients of a given polynomial can be extracted with 'coeff'. This function do not apply any simplification procedure (in contrast to 'apply', 'factor', and such):

The degree of a polynomial is extracted with 'hipow':

The following command returns a list containing the coefficients of the polynomial poly1:

Thus, the coefficients of powers of x in poly2 are:

and those of y:

With 'nroots' we can know the number of real roots a polynomial posseses in a given interval:

In order to try to find those roots exactly, we can use 'solve':

In many cases, we will not be able to compute roots in an exact way. With 'realroots' we obtain rational approximations:

Of course, in case the roots are integer or rational numbers, 'realroots' will give an exact answer:

Multiplicities are stored in the variable 'multiplicities', as it should be:

To compute real and complex roots in full generality, we can use 'allroots' (remark: when there are multiple roots, 'allroots' is not the best option. In those cases, it is better to try something like allroots(%i*p), if the polynomial has real coefficients):

Polynomial division using Euclid's algorithm is done with the aid of 'divide'. The output is a list of the form [quotient, remainder], whose elements can be extracted individually:

Maxima can also compute the greatest common divisor:

and the least common multiple:

Here, we express (x+a+1)^4 as a polynomial in x by grouping powers in x:

Now, let's consider some elementary single-variable calculus operations. For instance, limits:

The indeifnite integral of x^2 sin(x)^2:

By differentiating the output, we should get the original integrand, albeit possibly in a different (but equivalent) form:

In order to see the equivalence of both expressions, we can resort to some simplifying commands that apply trigonometric relations. In the first place, we have 'trigsimp', that will try to write the given expression with the minimal number of trigonometric functions:

THis result is expressed in terms of the cosine, and we want the sine. We then use 'trigexpand' to convert the cos(2x) term into factors containing sin(x) and cos(x):

and, finally, we simplify again (Maxima prefers the sine over the cosine for this, hence the trick):

Definite integration follows a similar syntax, adding the integration lower and upper limits, which can be infinite:

Numerical integration is also possible. Here we give the numerical value of int(sin(sinx)) between 0 and pi, using Romberg's method:

The 'diff' command admits several arguments to differentiate functions with any number of variables an arbitrary number of times. The next command also illustrates the fact that the apostrophe ' delays the evaluation of the function that it precedes:

The following command computes the Taylor expansion of the sin(x) function around x=0, up to degree 6:

Maxima is able to compute the Taylor development of a purely symbolic function. However, as the values of f are unknown, its derivatives are also unknown and must be declared in some way. A possibility is to use 'gradef', but a much better solution is to consider the package 'pdiff', which establishes a positional notation for the derivatives: given f=f(x), the n-th derivative with respect to x will be denoted f_(n)(x). If f=f(x,y) is a fucntion of two variables, f_(2,1)(x,y) will denote a third derivative of f, twice with respect to x and once with respect to y, etc. Thus, consider

The output of the command 'taylor' is a series, and as such it has properties which are different from those of a polynomial. If we want to use the 6th degree polynomial above, we need to make a truncation:

Finally, multiple integration works as one would expect:

 3 Elementary Number Theory

Maxima knows about the Golden Number (or Divine Proportion, Golden Constant, etc.)

However, Maxima does not establish its algebraic properties by default. We can do it by ourselves, declaring %phi as the only positive root of the polynomial x^2-x-1. The command 'tellrat' adds a number to the ring of algebraic integers already known by Maxima:

Here is the continued fractions expansion of %phi (of length 3):

Notice the length of the period. Things look better if we use a rational form for %phi:

Fibonacci numbers can be expressed in terms of %phi, a fact also implemented in stock Maxima:

It is also useful to have available commands for dealing with prime numbers:

We can determine whether a number is prime or not by applying some primality test. For integres greater than 341550071728321, probabilistic algorithms (such as Miller-Rabin and Lucas) are used:

The number of pseudo-primality tests to be run is given by the value of the variable 'primep_number_of_tests', here is its default value:

Maxima uses exact arithmetic when dealing with integers. The function 'factor', for instance, gives the prime factors decomposition of an integer:

It is possible to get a list with the prime factors and their exponents for later use:

It is known that the number ϕ(n)/n (where ϕ is Euler's totient function) approximates the probability that a positive integer chosen at random, be coprime with p1 , . . . , pk , where n = p1^e1 · · · pk^ek is the prime factorization of n. One could write first a program that selects m random numbers between 1 and 100m. The 'block' environment offers the possibility of defining some local variables (in this case, s and N) whose previous values will be stored and then erased, after exiting the block, they will be restored. Then, we can write a sequence of instructions (separated by commas) using the usual programming structures for loops, stop conditions, etc. The last instruction determines the output. In this regard, the following very simple function 'numgen' clearly illustrates these ideas:

For a list of such numbers, we can complete the fraction of those coprime with n (notice the appearance of the function 'gcd' for this):

Also notice the use of lambda functions here. Finally, we can check that P (m, n), for large values of m, approaches ϕ(n)/n:

 4 Graphics

There are two options available to create graphics in Maxima: the set of functions 'plot2d', 'plot3d' and the more powerful (and complex) 'draw2d', 'draw3d'. In each case, there exists a wrapper for wxMaxima that can be called just adding the prefix 'wx' to each command, resulting in graphics embedded into the wxMaxima worksheet (if only, say, plot3d is used, the graphics will appear in a separate window). Maxima uses the external grapher Gnuplot for these tasks, so it must be previously installed in your system.

Here is a simple drawing: the graph of the function  sin(exp(x)) with x running in [0,pi]:

Maxima tries to define a range of values for the y variable which are appropriate for the function at hand, but sometimes we need to make some adjustments. Consider, for example, the case of the function 1/(x+1)(x-1)^6:

The default range reflects the existence of a pole at x=-1 of order 1. However, the function has another pole at x=1 (of order 6) which can be made apparent by writing the function in an equivalent form:

It is clear that no single range of values will be adequate to accurately represent two poles of such disparate orders. A little bit of trial and error applied to [y,ymin,ymax] leads to a more or less satisfying compromise:

It is possible to define and draw piecewise-defined functions:

If we do not want the vertical lines that Gnuplot introduces in order to 'complete' the graph of the function, we can employ a little trick:

Here is a representation of the complex roots of z^n=1 in the plane:

Notice the 'same_xy' option here! It tells Maxima to use the same scale for the x and y axes: Maxima uses a 600x400 default canvas, so the original scales are in a 3:2 proportion, which would lead to a distorted graphics.

We have already remarked that Maxima uses an external grapher, Gnuplot, which posseses many options to adjust the graphical output. Here are some of them:

Gnuplot supports many types of lines. Their appearance is controlled by the options 'linetype' and 'dashtype', which are called using 'gnuplot_preamble', an instruction that allows us to use the Gnuplot syntax directly inside Maxima:

The contour diagram for the function sin(x+sin(y)). The option 'set cntrparam levels' determines the number of level sets to be shown:

This is better seen with a color code. Lighter areas belong to lower level sets:

The z axis marks accumulate creating a disappointing spot (the graphics is obtained by projecting the surface from above onto the plane z=0). A possibility for obtaining a clearer drawing consists in using 'draw3d' (previously loading the package 'draw'):

Here is a  3D representation of the same function (all of the following examples are better visualized if the instructions 'plot3d' and 'draw3d' are used, the reader can then click on the external Gnuplot window moving the mouse to rotate the graphics):

The 'plot3d' command works good with explicit surfaces of the form z=f(x,y), but if we need to create graphics for parametric curves or surfaces, it is better to use the several variants of 'draw' (in fact, 'draw' supports explicit, implicit as well as parametric drawings):

As an example, let us see how to represent a circumpherence explicitly, in parametric form and implicitly. We will take the opportunity to show other options, such as line color, line width, titles, etc.

Notice the use of 'proportional_axes=xy'! By default, draw uses a 600x400 canvas, so the plot of a circumpherence will appear distorted if we do not use that option (which has the effect of using the same scale in both axes, not the 3/2 proportion). It is the analog of the 'same_xy' option that we used in plot2d.

Here, we can see that the number of default 'ticks' used to construct the parametric curve is not enough. Maxima uses a default of 30 subdivisions of the interval, in this case [0,2pi], which gives the impression of a polygon, rather than a smooth curve. We can fix this by increasing the number of subdivisions with 'nticks':

Finally, an implicit drawing:

Implicit drawings require a lot of computational work, and should be avoided whenever possible. We can increase the quality of the output by means of 'ip_grid' (that will also increase the timing!)

Next, we generate a parametric surface, with the points x,y,z on the surface gievn as functions of the parameters u,v. It is important to notice the use of the 'wxplot_size' variable to match the size of the Gnuplot drawing with the size of the wxMaxima canvas in which we will draw it:

Here is a more complex parametric surface:

It is possible to combine several plots in one scene:

The 'draw' package can call a little library called 'worldmap' which contains cartographical data:

The next commands generate a two-scene graphics (better seen in an independent window, just replace 'wxdraw' below by 'draw'):

Let's see how to represent several 3D graphics together:

With 'draw' we can also create advanced 2D graphics. For instance, we can shadow the area between curves with 'filled_func'. Many more options can be found in the documentation.

 5 Data fitting

We have already had occasion to work with lists before. The following command creates a list containing the first 75 factorials (we chose not to show it with the symbol $ at the end of  the command):

Next, we take the logarithm of each element in the list and numerically evaluate the result. Most of Maxima built-in functions, like 'log' have the property of being 'listable', so they are applied individually to each element:

Individual elements of the list a:[a_1,...,a_n] can be extracted by specifying their indices, as in a[i]. Thus,

Let's check that the first number Lista[4] is the exponential of the second lista[4] (for this, we must express the floating numbers returned by the exponential as rationals.  This behavior is controlled by the environment variable 'bftorat'; the other variable 'ratprint', when set to 'false', avoids that some message telling that such or such number was converted to rational be emitted):

Now we generate a graphics representing the entries in lista against those in Lista. The fact that the points seem to distribute asymptotically near a straight line is a consequence of a mathematical result known as Stirling's formula: log(n!) behaves as nlogn, for big values of n.

Maxima has a function to fit data using the least squares technique, not surprisingly contained in the package 'lsquares'.

As Stirling's formula is asymptotic, we discard the first 24 values (which would lead to a distorted fitting)

First, our list of points must be converted into a matrix:

Then, the function "lsquares_estimates" fits the data in that matrix (whose first column is interpreted as containing the 'x' coordinates, while the second contains the 'y' ones) to a straight line y=mx+b, and gives us the values of the coefficients m,b:

Now we can plot the data and the fitting curve in the same graphics:

 6 Linear Algebra

Let us start by defining a matrix to play with:

Here we compute a basis of the null subspace (the kernel) of the linear mapping determined by C:

The result is given in the form of column vectors (neither matrices nor lists). Hence, the procedure to extract the basis vectors of Ker(C) is a little bit elaborated:

The dimension of ker(C) is directly given by 'nullity;:

Of course, we can also compute a basis of the column subspace (image subspace, Im(C))

and the corresponding dimension:

In the following command, we compute a basis of the orthogonal complement to Im(C):

We can check that these vectors and those in the basis B2 (of columnspace(C)) are orthogonal. Notice here how to define a matrix from a pre-existing function of two indexes as arguments:

Continuing with this idea, next we generate a 3x3 matrix whose (i,j) coefficient is 1/(i+j+1):

The inverse of M is

Multiplying the inverse by the original matrix we get the identity matrix (the usual . , not * , is used in Maxima to denote the matrix product):

The following command outputs a new matrix, with a modified diagonal. It is important to pay attention to the way the 3x3 identity matrix is written:

The determinant of this new matrix is precisely the characteristic polynomial of the original matrix:

The same computation can be done directly with the function 'charpoly':

With 'eigenvalues', we can find the numerical values of the eigenvalues of M (what a surprise!). Their multiplicities are also given inside the output, as a separate list:

The following commands generate a (pseudo) random matrix 1000x1000 (with real coefficients normally distributed around the mean 0 with variance 1):

In order to compute the eigenvalues of such a matrix, highly specialized routines, which can be found in the package 'lapack', must be used. In this case, we will need the command 'dgeev'.

Let us represent the output in the complex plane, to visualize the distribution of the eigenvalues:

The result is a numerical proof of Girko's Circular Law:

Symbolic matrices are also valid. In this example, we compute the eigenvalues and eigenvectors of a matrix of indeterminate coefficients:

Some particular matrices are already built-in, like Vandermonde matrices:

The determinant of such a matrix is computed by means of the classical formula:

Maxima can compute the Jordan normal form of a matrix directly. The related functions are contained in the package 'diag':

In this case, we have a diagonalizable matrix:

The output of the command 'jordan' is a list of eigenvalues with their multiplicities (or, in general, Jordan blocks). To display the matrix proper, we do

Let us consider the general case. Again, notice that the output of the functions is given in the form of lists:

Functions of matrices can be computed with the aid of the command 'mat_function', also contained in 'diag':

The Gram-Schmidt orthonormalization process with respect to the Euclidean canonical product in R^n (stored in the variable 'innerproduct') is carried on with the command 'gramschmidt', contained in the package 'eigen':

Let us check that the resulting vectors form an orthogonal system:

It is possible to define the scalar product to be used, and to pass it as an argument to 'gramschmidt'. In the following commands, we use this idea to construct the Legendre polynomials normalized so as to make their leading coefficient equal to 1:

A quick check:

 7 Differential Equations

Maxima can represent the slope field (isoclines) of a differential equation without problems. Here, we use the equation dy/dx=2x+y with x running in [-2,1] and y in [-4,6]:

The same command 'drawdf' (or its wxMaxima version) allows us to study the phase diagram of a first-order system (equivalently, a second-order differential equation), in this example
du/dt=v^2
du/dt=u^3

Problems involving initial conditions can be solved, too. The problem dy/dx=x+y with y(0)=2 uses the commands 'ode2' (for solving the equation) and 'ic1' (for applying the initial conditions):

There is a similar command, 'ic2', to deal with second-order equations:

Boundary-value problems for second-order differential equations are the target of 'bc2'. For instance, d²y/dx²+y=0 with conditions y(0)=0,y(pi)=1, has no solution:

The same equation with conditions y(0)=0,y(pi/2)=1 has a unique solution:

Here is the phase diagram of a pendulum with the trajectories through some points selected with the function 'solns_at':

Ordinary differential equations with parameters are admissible, too. Here we compute the solution of y''(x)-ky(x)=1, assuming k>0:

Had we not used the assumption on k above, Maxima would ask us about the sign of k. The form of the solution depends crucially on the sign of k, as we can check:

To numerically solve a second-order differential equation, Maxima offers the function 'rk', which implements the Runge-Kutta method of order four. However, this function only accepts first-order differential equations or systems, so in order to solve an equation like y''(x)+sin(x)^2 y'(x)+y(x)=cos(x)^2, with initial conditions y(0)=1, y'(0)=0, we must write it as a first-order system first. The procedure is standard: introduce an auxiliary function u(x) and write
y'(x)=u(x)
u'(x)=-sin(x)^2 u(x)-y(x)+cos(x)^2
Now we can apply 'rk'. For instance, to compute y(x) with x running in [0,20]:

The result is a list of sublists having the form [xi,yi,ui]. We are interested in the graphical representation of the solution y(x), that is, the first two elements of each sublist, that we can extract as follows:

 8 Bessel Functions

Bessel functions are implemented in stock Maxima, but they are treated symbolically (that is, they are not evaluated by default) except by certain functions like 'specint'. This is the same that happens when you write

Probably you would expect to see a^3+3ba^2+3b^2a+b^3, but this is not the output, because Maxima treats symbolic expressions in a symbolic manner, so it sees (a+b)^3 as a whole. If you want the expansion of this expression, you have to ask for it explicitly:

Well, as mentioned, the same happens in the case of Bessel functions:

If we want a numerical result, we must ask for it:

To draw the Bessel functions, we can use the evluation operator on functions, given by '' (two simple apostrophes). This forces the function to evaluate the result, not just to store it. Now everything is easy:

The graphics appearing in every textbook:

We could look for roots in an interval [a,b] with the command 'find_root' but we must keep on mind that the computation is stopped once the first root is found, so it is not possible to apply it to list all the roots in a big interval such as [0,10]. For instance, in the case of J0 we could do something like

The calculations can be simplified if we recall the Sturm theory, which says that the distance between two consecutive zeros tends to pi. Thus, we can look for zeros in an interval by dividing it into segments of length close to pi. In the case of the zeros of J0 over the interval [0,10] we would do

A simple function for computing the zeros of J0 in [a,b] can be written along these lines:

so the above example would read now

Our little function can answer questions such as: what are the first 40 zeros of J0? The result is given in the form of a list, named 'first40' so we can use it later:

Of course, we can numerically prove the Sturm prediction stating that the interval between consecutive roots tends to pi:

As an easy exercise, the reader will have no problem in adapting the preceding computations to the case of other Bessel functions J_nu and Y_nu.

 9 Animations

Maxima, when used in conjunction with its interface wxMaxima, allows the user to create very basic animations depending on a parameter (the 'slider'). Next, we show an illustration of the tangent line to a function at a point (to see the animation, just click on the graphics with the left button of the mouse, or with the right one and select "Start Animation"):

Here we animate Archimedean spiral %rho=%theta with %theta running from 8%pi to 10%pi:

Let us consider some more advanced animations. The transformation of a rectangular strip into an annulus:

and a transformation of the annulus into a Moebius band:

An animation illustrating the definition of the circular (trigonometric) functions:

and the torus as a surface of revolution:

As our final example, here is an animation of the motion of a circular membrane (the expression below was obtained by solving a Dirichlet problem for the Laplacian on the disk, which involves the Bessel functions):

 10 Lagrangian and Hamiltonian of the Electromagnetic Field

As an application to Physics, in this section we study the construction of the Lagrangian for the classical electromagnetic field and, from it, the canonical Hamiltonian. The corresponding equations of motion are nothing but the Lorentz force equations.

Let's fix the space dimension

Define the generalized coordinates

Now we declare the time dependency

Define the velocities:

Velocities are the derivatives of generalized coordinates:

The potentials (depending on the generalized coordinates):

Let's check what do we have up to now:

Now comes the Lagrangian:

From L we can define the conjugate momenta:

For consistency, we can also store the momenta in a list of variables [p1,p2,p3]. Notice the use of the macro operator '::' to make the assignment and the evaluation of the left-hand side all at once:

Notice that

but these are actually the momenta: we only need to evaluate to see it explicitly

Later on, we will see the convenience of defining the list of variables p=[p1,p2,p3]. Now, we use the momenta to get the Euler-Lagrange equations:

As a quick check, we can see that these equations are actually those corresponding to the Lorenz force m dv/dt =ev x rotA + eE. First, we define our own version of the curl (or 'rotational', hence the name), taking benefit of the particular dimension n=3:

This works for vector fields whose dependence has been declared by 'depends':

We will need the cross product operator '~', which is contained in the package 'vect':

Equations for the Lorenz force:

Those are exactly the Euler-Lagrange equations obtained before:

Let us see if this Lagrangian is regular. No need to reinvent the wheel here: the package 'linearalgebra' already contains a command 'Hessian'

As m>0, we see that the Lagrangian of the electromagnetic field is regular. We can proceed to find the Hamiltonian through the Legendre transform. We will use the momenta pj as new variables, using the list p=[p1,p2,p3] defined earlier. First of all, we need to express the velocities in terms of the momenta

Below we define the Hamiltonian, and simplify it to get the standard form

In order to get the standard form of the Hamiltonian, we need to complete squares. We define a little function for that:

 11 Signal Analysis

Some basic signal analysis can be done in Maxima. We will need to generate random numbers:

Number of samples-1:

Sampling points (2^10 of them, evenly spaced):

Let us add some white noise to a sinusoidal wave (of frequency w=211):

The result seems to be a completely random signal:

However, a frequency analysis (using the fast Fourier transform) reveals that inside the signal there is a regular, periodic wave of frequency... about w=211!

Now, we superimpose three rapidly oscillating signals of frequencies w0=1920, 2w0 and 4w0:

And repeat the procedure:

The three peaks are identified, with their relative separation, but their locations are displaced. This is due to the fact that the list u only contains 513 values (1024/2+1). To put the frequencies in their absolute abcissas, we must sample above the highest frequency: