iconEuler Home

Symbolic Expressions and Maxima

This is a short introduction to Maxima for Euler. For more information
on Maxima view their homepage at Sourceforge, or the reference
delivered and installed with Euler.

http://maxima.sourceforge.net/
See: ../reference/maximacore
See: 00 - Overview of the Syntax of Euler

The preferred method to communicate with Maxima are symbolic
expressions and functions. Symbolic functions start with "&" and can
be used at any place in Euler, where expressions are allowed.
>plot2d(&diff(x^4+exp(-x^2),x),-1,1):

10 - Introduction to Maxima

They can also be used directly, and print to the notebook in Maxima's
formula display.
>&diff(x^4+exp(-x^2),x)
                                          2
                              3        - x
                           4 x  - 2 x E

It is possible to output a symbolic term with Latex, if you have Latex
installed. However this will slow down the execution of your
notebooks.

It might be better to use "maxima: formula" in a comment. Here is an
example.

maxima: 'integrate(x^2*log(x),x)=factor(integrate(x^2*log(x),x))+c

The code for this was

 maxima: 'integrate(x^2*log(x),x)=factor(integrate(x^2*log(x),x))+c

Here is a direct output of a symbolic expression in the usual Euler
output.
>$factor(diff(x^4*exp(-x^2),x))

10 - Introduction to Maxima

Symbolic expressions can be stored in variables accessible for other
symbolic expressions.
>fx &= log(x)/x
                                log(x)
                                ------
                                  x

It helps to be aware that fx is defined for Euler in a string, which
prints as a symbolic expression.
>fx
                                log(x)
                                ------
                                  x

To print it as a string, you can use string concetanation.
>""|fx
log(x)/x
To get a Latex representation, there is the tex command.
>tex(fx)
\frac{\log x}{x}
The Latex output can be produced with $fx.
>$fx

10 - Introduction to Maxima

Now we can plot the function in Euler as usual. This works, since the
function is just an expression in x, stored in a string.
>plot2d(fx,1,10):

10 - Introduction to Maxima

Symbolic expressions are evaluated in Maxima, and the result is just
an Euler string. We use this result string in the function solve to
find the numerical value of the maximal point close to the initial
guess 3.

Note that in this command "solve" is the Euler function "solve".
>solve(&diff(fx,x),3)
2.71828182846
For this simple function fx, we can also use Maxima to solve the
value exactly.

In the following command, we use the Maxima command "solve".
>&solve(diff(fx,x))
                               [x = E]

We can apply other numerical methods to the function.
>romberg(&fx,1,E)
0.5
This is indeed the correct integral. Maxima can compute the integral
for this simple function.
>&integrate(fx,x,1,E)
                                  1
                                  -
                                  2

We can easily define an expression for the second derivative. The
right hand side of &= is parsed through Maxima before the definition.
>d2fx &= factor(diff(fx,x,2))
                             2 log(x) - 3
                             ------------
                                   3
                                  x

Find the the turning point numerically.
>solve(d2fx,E,4)
4.48168907034
And exactly, which is possible in this case.
>&solve(d2fx,x)
                                    3/2
                              [x = E   ]

We can also define a symbolic function in Euler and Maxima with the
same expression. The function definition uses &= in this case. 

The right hand side is evaluated in Maxima before the definition of
the function. In this case, it simply calls the symbolic expression
fx.
>function f(x) &= fx
                                log(x)
                                ------
                                  x

Maxima can use this function, as well as Euler.
>&f(4), f(4),
                                log(4)
                                ------
                                  4

0.34657359028
Here is a function for the tangent line of f in x0.
>function T(x,x0) &= factor(f(x0)+(x-x0)*at(diff(f(x),x),x=x0))
                  2 x0 log(x0) - x log(x0) - x0 + x
                  ---------------------------------
                                   2
                                 x0

Note that &= parses the function body through Maxima, before it
defines the function.
>type T
function T (x, x0)
useglobal; return (2*x0*log(x0)-x*log(x0)-x0+x)/x0^2 
endfunction
Let us plot the function and its tangent in one plot.
>plot2d("f(x)",E,10); plot2d("T(x,E^(4/2))",add=1,color=4):

10 - Introduction to Maxima

If we try to find a tangent, which goes through (0,-2), we fail in
Maxima. There is no closed solution for this.
>&solve(T(0,x)=-2,x)
                               1 - 2 log(x)
                          [x = ------------]
                                    2

However, in Euler we can easily find a numerical solution.
>solve(&T(0,x)+2,1)
0.766248608162

Direct Input for Maxima

Another way to use Maxima in Euler is to send commands directly to
the Maxima system, which runs as a process and communicates with
Euler. Prepend "::" or ":" in front of the command to send it to
Maxima.

":: " is the compatibility mode
": " sends the command verbatim to Maxima

The compatibility mode makes the Maxima syntax similar to the Euler
syntax. E.g., it uses ":=" to set variables instead of ":", and uses
";" to seperate commands and not print results, and "," to print
results, just as Euler does. The direct mode is for experts in
Maxima.
>:: 30!
                  265252859812191058636308480000000

>:: factor(%) // Factor the previous result
                   26  14  7  4   2   2
                  2   3   5  7  11  13  17 19 23 29

With "%", it is possible to refer to the previous output, just as in
Euler. However, to access results over several lines it is better to
use variables.

As in Euler, "//" starts a comment.

A line can contain several Maxima commands. The commands are sent to
Maxima one after the other. Separate the commands with "," or with
";".
>:: expand((1+x)^10), factor(diff(%,x))
        10       9       8        7        6        5        4
       x   + 10 x  + 45 x  + 120 x  + 210 x  + 252 x  + 210 x
                                                  3       2
                                           + 120 x  + 45 x  + 10 x + 1


                                       9
                             10 (x + 1)

In case, you wish to use the direct mode, remember that "$"
suppresses the output of a command, ":" is used to set variables, and
"," is used for flags to a command. The delimeter ";" prints the
output of a command.

By default the direct mode is disabled. So you will probably get an
error on your Euler system.
>: expr:1+x$ (1+x)^10; %, expand
                                     10
                              (x + 1)


         10       9       8        7        6        5        4
        x   + 10 x  + 45 x  + 120 x  + 210 x  + 252 x  + 210 x
                                                  3       2
                                           + 120 x  + 45 x  + 10 x + 1

Maxima Mode

It is also possible to switch to Maxima mode permanently. All
commands are Maxima commands from that point on. Euler can be
configured to start in Maxima mode automatically.

We use this mode to demonstrate some features of Maxima.

To use Euler commands in Maxima mode, prepend the command with "euler
...".
>maximamode on
Maxima mode is on (compatibility mode)
Note that variables are set with := in Maxima. Here we set a=3,
then evaluate an expression in a, then kill the variable a.
>a:=3; (1+a)^3, remvalue(a);
                                  64

Another way to evaluate an expression with a value is to locally set
the variable. The flag "with ..." does this. Append it to the
command.
>(1+a)^3 with a=3
                                  64

The third method to evaluate an expression with a value is to
define a Maxima function.
>function f(a) := (1+a)^3
                                          3
                           f(a) := (1 + a)

>f(3)
                                  64

Simplification

Maxima does not automatically simplify completely. You need to
say, what you want. In the following example, we go back and forth
from one form to another of a rational expression.
>2/a + a + 1/a^2, factor(%), expand(%)
                                  2   1
                              a + - + --
                                  a    2
                                      a


                              3
                             a  + 2 a + 1
                             ------------
                                   2
                                  a


                                  2   1
                              a + - + --
                                  a    2
                                      a

Another way to simplify a rational function is to apply the command
ratsimp to the evalutation.

Such flags have to be appended using | ...
>2/a + a + 1/a^2 | ratsimp
                              3
                             a  + 2 a + 1
                             ------------
                                   2
                                  a

Appending the command after a "|" is equivalent to a direct command
for many simplifications.
>ratsimp(2/a+a+1/a^2)
                              3
                             a  + 2 a + 1
                             ------------
                                   2
                                  a

"ratsimp" and "factor" will also cancel common factors in the
numerator and denominator.
>(a^3+3*a+3*a^2+1)/(a+1), factor(%)
                          3      2
                         a  + 3 a  + 3 a + 1
                         -------------------
                                a + 1


                                      2
                               (a + 1)

Some other commands can simplify more complex expressions. You
might want to try "radcan" and "trigsimp".
>radcan(log(100)/log(10))
                                  2

>trigsimp(sin(x)^2+cos(x)^2)
                                  1

"simpsum" is a command, which simplifies and evaluates a sum.
>sum(k^2,k,1,n) | simpsum | factor
                         n (n + 1) (2 n + 1)
                         -------------------
                                  6

Calculus

Of course, Maxima can differentiate and integrate too.
>diff(x^x,x)
                            x
                           x  (log(x) + 1)

Here is the second derivative.
>diff(x^x,x,2)
                       x             2    x - 1
                      x  (log(x) + 1)  + x

You can use "at" to evaluate the derivative at a specific point.
>at(diff(x^x,x,2),x=1)
                                  2

However, it might be easier to use the funciton "diffat", which
is an extension of Maxima in Euler.
>diffat(x^x,x=1,2)
                                  2

Other extensions are gradient, hessian and jacobian. These functions
return Maxima matrices. See below for more information.
>gradient(x^y,[x,y])
                          y - 1     y
                        [x      y, x  log(x)]

>hessian(x^2+y^2+x*y,[x,y])
                               [ 2  1 ]
                               [      ]
                               [ 1  2 ]

>jacobian([x^y,y^x],[x,y])
                       [  y - 1      y        ]
                       [ x      y   x  log(x) ]
                       [                      ]
                       [  x            x - 1  ]
                       [ y  log(y)  x y       ]

If you are in a hurry, omit the variable list, using grad(f) or
hesse(f).
>grad(x^2+y^2)
                              [2 x, 2 y]

For the following, we switch off Maxima mode.
>maximamode off
Maxima mode is off

Calculus with Symbolic Expressions

Since strings, which contain expressions in x can be evaluated at
specific values of x, this makes the evaulation of a Maxima result
easy.
>&diff(x^x,x,2)(4)
1521.76659915
Alternatively, the string can be evaluated with global values.
>x:=4; &diff(x^x,x,2)()
1521.76659915
There is a shortcut for this. &:expression evaluates the expression
and uses the result directly as an Euler command.
>x:=4; &:diff(x^x,x,2)
1521.76659915
Alternatively, we can store the expression into variables of Euler
and Maxima, and evaluate later.
>x:=4; expr &= diff(x^x,x,2); expr()
1521.76659915
The symbolic gradient function returns a symbolic vector.
>&gradient(x^2-y^2,[x,y])
                             [2 x, - 2 y]

The numerical value of the Hessian matrix is computed by mxmhessian.
We evaluate the result at x=y=E.
>&hessian(x^y-y^x,[x,y])(E,E)
     -5.57494152476                   0 
                  0       5.57494152476 
The same method works for integrals.
>&integrate(x*sin(x),x,0,1)()
0.30116867894

Tables with Maxima

An interesting command is mxmtable, which creates a table of values
using Maxima. It prints the values and a plot of the values. The
second argument is the variable name.
>mxmtable("sum(1/k,k,1,n)","n",1,10);
                  1 
                1.5 
      1.83333333333 
      2.08333333333 
      2.28333333333 
               2.45 
      2.59285714286 
      2.71785714286 
      2.82896825397 
      2.92896825397 

10 - Introduction to Maxima

Even a table of two variables can be created.
>mxmtable(["2^-n","3^-n"],"n",0,10,print=0);

10 - Introduction to Maxima

Variables in Euler and Maxima

Note that all variables in Maxima are different from the variables in
Euler. But symbolic expressions are known to both worlds.
>expr &= x/(1+x^2)
                                  x
                                ------
                                 2
                                x  + 1

For numerical values and small matrices, we can use ::= to set a
value for Maxima and Euler. This value can be used in all symbolic
expressions.
>v &:= 1:10; &v[10]
                                  10

For large matrices, we need to use mxmset. 

The subtle difference is that v will always be a Maxima matrix now,
not just a Maxima vector. This there must be two indices, for row and
column.
>mxmset("v",1:1000); &v[1,1000]
                                 1000

Newton Method via Maxima

Symbolic algebra is useful to help numerical algorithms in Euler.
We can use the Newton method, and let Maxima compute the derivative
of the function for us.
>longestformat; newton("x^x-2",&diff(x^x-2,x),2)
      1.559610469462369 
The utility function "mxmnewton" makes this easier.
>mxmnewton("x^x-2",2)
      1.559610469462369 
To check, we use the bisection method.
>bisect("x^x-2",1,2)
      1.559610469462314 
The interval Newton method provides a guaranteed inclusion of
the solution. Note that the bisection does not have this accuracy.
>mxminewton("x^x-2",~1,2~)
    ~1.5596104694623687,1.55961046946237~ 
Maxima can also be used for the Newton method in several dimensions.

The utility function "mxmnewtonfxy" uses the Newton method to
find a point where both expressions are 0.
>longformat; mxmnewtonfxy("x^2+y^2-10","x+y-4",[1,2])
[ 1  3 ]
Let us demonstrate how to do this without utility functions.

First we define the function as a symbolic expression.
>fexpr &= [x^2+y^2-10,x+y-4]
                        2    2
                      [y  + x  - 10, y + x - 4]

Then we define the function, using the vector syntax for the function
arguments. In general, we need a function, which takes a vector and
returns a vector. We define it as a symbolic function.
>function f([x,y]) &= fexpr
                        2    2
                      [y  + x  - 10, y + x - 4]

Moreover, we need a function, which takes a vector, and returns the
Jacobian. To compute the derivatives, we use Maxima at compile time.
>function Df([x,y]) &= jacobian(fexpr,[x,y])
                             [ 2 x  2 y ]
                             [          ]
                             [  1    1  ]

Indeed, the function is defined as expected.
>type Df
function Df ([x, y])
useglobal; return matrix([2*x,2*y],[1,1]) 
endfunction
Again, it is easier to use symbolic expressions.
>function Df([x,y]) &= jacobian(fexpr,[x,y])
                             [ 2 x  2 y ]
                             [          ]
                             [  1    1  ]

Now we can start the Newton method.
>newton2("f","Df",[1,4])
[ 1  3 ]
Or we can start an interval Newton method, which provides a guaranteed
inclusion.
>inewton2("f","Df",[1,4])
[ ~0.999999999999998,1.000000000000003~
~2.999999999999997,3.000000000000003~ ]

Maxima Questions

Sometimes, Maxima is asking questions. If this happens, enter
the answer in the next input line.
>:: integrate(x^n,x)
Is  n + 1  zero or nonzero?

Here, we choose nonzero.
>:: nonzero
                                 n + 1
                                x
                                ------
                                n + 1

The mxm command does not answer questions. You need to use the
assume feature of Maxima to setup everything.
>&assume(n>0); &integrate(x^n,x)
                                 n + 1
                                x
                                ------
                                n + 1

Eigenvalues

Let us generate a matrix in Euler and Maxima.
>shortformat; A:=2*id(3); A:=setdiag(A,-1,1); A::=setdiag(A,1,1)
Syntax error in expression, or unfinished expression!

Error in :
shortformat; A:=2*id(3); A:=setdiag(A,-1,1); A::=setdiag(A,1,1)
                                               ^
Compute the eigenvalues in Maxima.
>ev &= eigenvalues(A)
              [[2 - sqrt(2), sqrt(2) + 2, 2], [1, 1, 1]]

This returns a vector with two elements, the eigenvalues and their
multiplicities. We extract and evaluate the eigenvalues.
>&:ev[1]
[ 0.58579  3.4142  2 ]
This is the same result as numerically in Euler.
>re(eigenvalues(A))
[ 0.58579  2  3.4142 ]

Integration

Here is an integral, where Maxima uses the function erf for the
result.
>&integrate(exp(-x^2/2),x,-1,1)
                                            1
                    sqrt(2) sqrt(pi) erf(-------)
                                         sqrt(2)

We can evaluate the last expression numerically.
>longestformat; &integrate(exp(-x^2/2),x,-1,1); %()
      1.711248783784292 
This can be done in one step.
>&:integrate(exp(-x^2/2),x,-1,1)
      1.711248783784292 
We can also check with the Romberg method integrating the Gauß
distribtion.
>romberg("exp(-x^2/2)",-1,1), longformat;
       1.71124878378431 
The line integral of the Gauß normal distribution is exact in
Maxima, of course.
>&integrate(exp(-x^2/2),x,minf,inf)/sqrt(2*pi)
                                  1

Partial Fractions

It is easy to get a partial fraction decomposition in Maxima.
>&partfrac((x-1)/(x^2*(x+1)),x)
                               2     2   1
                           - ----- + - - --
                             x + 1   x    2
                                         x

However, let us try to do that by hand.
>&remvalue(A,B,C); &ratsimp(A/x+B/x^2+C/(x+1)), &expand(num(%))
                     2                   2
                    x  C + (x + 1) B + (x  + x) A
                    -----------------------------
                                3    2
                               x  + x


                      2                2
                     x  C + x B + B + x  A + x A

Now we compare coefficients and solve for A, B and C.
>&solve([coeff(%,x,0)=-1,coeff(%,x,1)=1,coeff(%,x,2)=0],[A,B,C])
                     [[A = 2, B = - 1, C = - 2]]

Then substitute the solutions back.
>&A/x+B/x^2+C/(x+1) with %[1]
                               2     2   1
                           - ----- + - - --
                             x + 1   x    2
                                         x

Curve Length

Here is a curve in the form of an "8".
>plot2d("sin(x)","sin(x)*cos(x)",xmin=-pi,xmax=pi,r=1):

10 - Introduction to Maxima

We compute the curve differential using Maxima.
>ds &= trigsimp(sqrt(diff(sin(x),x)^2+diff(sin(x)*cos(x),x)^2))
                             4           2
                   sqrt(4 cos (x) - 3 cos (x) + 1)

Now we compute the length of the curve numerically with Euler.
>longestformat; romberg(ds,-pi,pi)
      6.097223470105731 

Maxima Packages

Maxima for Euler can load packages.
>&batch(solve_rec);
We solve the Fibonacci recursion numerically.
>sol &= solve_rec(h[n+2]=h[n+1]+h[n],h[n])
                           n          n                n
              (sqrt(5) - 1)  %k  (- 1)    (sqrt(5) + 1)  %k
                               1                           2
         h  = ------------------------- + ------------------
          n               n                        n
                         2                        2

To get h1=1, h2=1, we solve for the constants.
>&solve([subst(0,n,rhs(sol))=0,subst(1,n,rhs(sol))=1],[%k[1],%k[2]])
                               1              1
                  [[%k  = - -------, %k  = -------]]
                      1     sqrt(5)    2   sqrt(5)

And substitute the constant back to the solution. The right
hand side is our desired formula.
>&sol with %[1]; sol &= rhs(%)
                             n                n      n
                (sqrt(5) + 1)    (sqrt(5) - 1)  (- 1)
                -------------- - ---------------------
                           n                   n
                  sqrt(5) 2           sqrt(5) 2

Now we can plot the logarithmic growth of the Fibonacci numbers.
>n:=1:100; m:=mxmeval("sol",n=1:100); plot2d(n,log(m)):

10 - Introduction to Maxima

Limits

The following limit is no problem.
>sx &= (1+1/x)^x, &limit((1+1/x)^x,x,inf)
                                1     x
                               (- + 1)
                                x


                                  E

If we want to estimate the convergence, we need to use Taylor series
methods. Maxima does that automatically.
>&limit(x*(sx-E),x,inf)
                                   E
                                 - -
                                   2

Even finer:
>&tlimit(x*(x*(sx-E)+E/2),x,inf)
                                 11 E
                                 ----
                                  24

Interpolation

Let us interpolate with a rational function.
>p &= (a*x+b)/(1+c*x)
                               a x + b
                               -------
                               c x + 1

To find a,b,c, we use the solve command of Maxima.
>&solve([(p with x=0)=1,(p with x=1)=-1,(p with x=2)=2],[a,b,c])
                             7               5
                     [[a = - -, b = 1, c = - -]]
                             6               6

Substitute back into p.
>fres &= ratsimp(p with %[1])
                               7 x - 6
                               -------
                               5 x - 6

We can also make a symbolic function out of this.
>function f(x) &= fres
                               7 x - 6
                               -------
                               5 x - 6

As you see in the following listing of f, the Maxima result has
been substituted into the function.
>type f
function f (x)
useglobal; return (7*x-6)/(5*x-6) 
endfunction
Test the interpolation values.
>longformat; f([0,1,2])
[ 1  -1  2 ]
If we need the Maxima result inside a string, we can use "&:...".
>function F(x) := romberg("&:fres",0,x)
>type F
function F (x)
useglobal; return romberg("(7*x-6)/(5*x-6)",0,x) 
endfunction
>F(0.5)
0.441281679648

LA-PACK

Since Maxima includes an implementation of LAPACK, we can use
that too..
>&load(lapack);
Let us try a random matrix.
>A:=normal(10,10);
First compute the eigenvalues in Euler.
>longestformat; sort(abs(eigenvalues(A)))
Column 1 to 2:
     0.1848629591742844      0.8529885918643476 
Column 3 to 4:
     0.8529885918643478       1.406747207126215 
Column 5 to 6:
      1.798254702671476       2.510395156517675 
Column 7 to 8:
      2.510395156517675       2.781614516683229 
Column 9 to 10:
      3.238847068210578       3.238847068210578 
Then in Maxima.
>mxmset(A,large=1); sort(abs(&:dgeev(A)[1]))
Column 1 to 2:
       0.18486295917428      0.8529885918643514 
Column 3 to 4:
     0.8529885918643514       1.406747207126214 
Column 5 to 6:
      1.798254702671493       2.510395156517692 
Column 7 to 8:
      2.510395156517692       2.781614516683248 
Column 9 to 10:
      3.238847068210579       3.238847068210579 
For badly conditioned matrices, Euler produces consistent results
with two methods.
>sort(jacobi(hilb(10)))
Column 1 to 2:
    0.01591324125174927       3.299660244880538 
Column 3 to 4:
      312.5985833767118       17889.84958322093 
Column 5 to 6:
      688491.8730021311       18741836.42660227 
Column 7 to 8:
      368416954.6185145       5202868202.111884 
Column 9 to 10:
      49919602009.29112       255023613679.9151 
>sort(abs(eigenvalues(hilb(10))))
Column 1 to 2:
    0.01591288576042408       3.299661583870189 
Column 3 to 4:
      312.5985816454515       17889.84958314523 
Column 5 to 6:
      688491.8730022628       18741836.42660274 
Column 7 to 8:
      368416954.6185171       5202868202.111885 
Column 9 to 10:
      49919602009.29113       255023613679.9151 
LAPACK produces the same result.
>mxmset("A",hilb(10),large=1); sort(abs(&:dgeev(A)[1]))
Column 1 to 2:
      0.015912780426876       3.299660422483491 
Column 3 to 4:
      312.5985841071862        17889.8495843014 
Column 5 to 6:
         688491.8730018       18741836.42660226 
Column 7 to 8:
      368416954.6185204       5202868202.111908 
Column 9 to 10:
      49919602009.29111        255023613679.915 
Here, the singular values agree with the eigenvalues.
>{u,v,w}:=svd(hilb(10)); sort(v),
Column 1 to 2:
    0.01591320516192639       3.299660128123471 
Column 3 to 4:
      312.5985834517168       17889.84958279682 
Column 5 to 6:
       688491.873001818        18741836.4266019 
Column 7 to 8:
      368416954.6185153        5202868202.11189 
Column 9 to 10:
       49919602009.2911       255023613679.9151 
In Maxima, we can use dgesvd only if we convert the Hilbert matrix
to float. The Hilbert matrix in Euler is scaled by hilbertfactor
to get an exact representation in floating point.
>sort(&:dgesvd(float(hilbert_matrix(10)),true,true)[1]*hilbertfactor)
Column 1 to 2:
    0.01591382668540789       3.299658163972339 
Column 3 to 4:
      312.5985839680668       17889.84958202317 
Column 5 to 6:
      688491.8730032905       18741836.42660133 
Column 7 to 8:
      368416954.6185104        5202868202.11185 
Column 9 to 10:
      49919602009.29123       255023613679.9151 

Euler Home