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):
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))
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
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):
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):
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
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
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
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
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
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
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
Even a table of two variables can be created.
>mxmtable(["2^-n","3^-n"],"n",0,10,print=0);
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
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~ ]
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
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 ]
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
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
Here is a curve in the form of an "8".
>plot2d("sin(x)","sin(x)*cos(x)",xmin=-pi,xmax=pi,r=1):
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 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)):
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
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
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