# Maxima Crash Course

## Maxima as a programming environment

Computer algebra systems generally follow a functional programming paradigm. A number of helper and higher-order functions are available in
Maxima. 
Map and apply are typical of the functional paradigm.

map("=",[x,y],[a,b]);

Example: Use "apply" to sum the elements of a list.

apply("+",[1,2,3]);

Makelist creates a list using a general form by running through an index.

makelist(j^2,j,1,10);

Evaluation vs substitution.

f: sin(x^5)$

ev(f,x=5);

subst(5,x,f);

subst(x,5,f);

subst(cos(t),5,f);

g: sqrt(x^2 + 2*x + 3);

ratsubst(14+y^2,x^2 + 2*x,g);

subst(14+y^2,x^2+2*x,g);

The last two show that subst only works on atoms or complete subexpressions.

atom(x);

atom(cot(x));

The quote operator prevents evaluation of an expression but not simplification.

'sin(0);

Maxima consider the above simplification rather than evaluation.

f(x) := 2*x;

'f(2);

Get documentation on a particular function.

? map;

## Calculus

Evaluate limits.

limit((1+1/n)^n,n,inf);

Derivatives.

f: sin(2*x) + exp(-x)$

diff(f,x);

Partial derivatives and the total differential. 

g: exp(-x*y)*cos(sqrt(y))$

diff(g,y);

diff(g);

One could use the total differential and 'solve' function to do implicit differentiation.

expr: x^2 + x*y^2 = 9$

diff(expr);

solve(%,del(y));

%/del(x);

Antiderivatives and definite integrals analytically.

f: cos(x)*sin(x)^2;

g: exp(-t)/t;

integrate(f,x);

integrate(f,x,0,1);

integrate(g,t);

This integral was recently reported as a bug case.

f: log(5 - 4*cos(x))$

integrate(f,x,0,%pi);

quad_qags(f,x,0,%pi);

first(%);

guess_exact_value(%);

The quote operator has a cousin which is frequently useful. Consider the following example of defining one function in terms of a derivative.

y(t) := t^2;

h(t) := diff(y(t),t);

h(2);

h(t) := ''(diff(y(t),t));

h(2);

This works because the quote-quote operator causes evaluation when the function is defined, so the derivative is computed before substituting t=2. 

How to calculate a path integral? Compute along γ from 0 to 1:

f: 4*x^3*y$

γ: [1-3*t,2-3*t]$

dγ: diff(γ,t);

dS: sqrt(dγ . dγ);

map("=",[x,y],γ);

sublis(%,f);

integrate(%*dS,t,0,1);

## Matrices and linear algebra

load("linearalgebra")$

Define matrices.

A: matrix([2,4],[1,2]);

B:  matrix([6/7,4,29],[11,0,1],[0,3,0]);

I: ident(5);

The dot operator " . " is used for matrix multiplication and the standard inner product.

A . A;

The ordinary multiplication operator performs the Hadamard product (element-wise).

D: identfor(B)*[x,1,1/y];

Difference between matrix^n and matrix^^n:

A^2;

A^^2;

There are two available functions to compute the determinant of a square matrix. "newdet" is supposed to be faster.

determinant(B);

Matrix inversion.

invert(B);

The Moore-Penrose pseudoinverse is available via the linearalgebra share package.

moore_penrose_pseudoinverse(A);

rank(A);

The matrix exponential.

matrixexp(A);

More generally, one may compute an analytic function of a square matrix, provided the corresponding power series of the matrix argument is
convergent. 

load("diag")$

mat_function(sin,A);

Eigensystem of a matrix.

load("eigen")$

eigenvalues(A);

eigenvectors(A);

Rather than have a separate function to diagonalize a matrix, Maxima just has the Jordan normal form via the diag package.

load("diag")$

jordan(B);

dispJordan(%);

There are some high-performance numerical linear algebra functions via LAPACK (Linear Algebra PAckage), which was translated to Common Lisp
from FORTRAN. (Matrices for these functions should have float or integer entries).

load("lapack")$

dgesvd(float(B));

Compare to analytically calculating the singular values of B. 

transpose(B) . B;

eigenvalues(%);

Should there be complex eigenvalues? These are eigenvalues of a real symmetric matrix, so they should have zero imaginary part. The reason this
occurs is that imaginary numbers may arise in using Cardano's formula for the cubic equation. We show that the imaginary part is indeed zero:

trigreduce(rectform(%));

float(%);

And the singular values are the square roots:

sqrt(%);

## Some fun: Lorenz attractor

bounds(x,y,z) := (
   makelist(
        [apply(min, l),apply(max, l)],l,[x,y,z]
    ) 
)$

round_bounds(L) := (
makelist(
    [floor(L[i][1]), ceiling(L[i][2])], i, 1, 3)
)$

kill(x,y,z)$
[sigma, rho, beta]: [10, 28, 8/3]$
eq: [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]$
sol: rk(eq, [x, y, z], [1, 0, 0], [t, 0, 50, 1/80])$
len: length(sol)$
x: makelist(sol[k][2], k, len)$
y: makelist(sol[k][3], k, len)$
z: makelist(sol[k][4], k, len)$
L: round_bounds(bounds(x,y,z))$
wxanimate_framerate:20$
wxplot_size:[840,700]$
with_slider_draw3d(
    k,makelist(i,i,1,len,20),
    title = concat("Lorenz Attractor"),
    points_joined=true,
    color = midnight_blue,
    point_type = -1, 
    xrange = [L[1][1],L[1][2]], yrange = [L[2][1],L[2][2]], zrange = [0,L[3][2]],
    points(firstn(x,k),firstn(y,k),firstn(z,k)),
    view = [75,35],
    axis_3d = true
)$