cadadr.org / maxima

Maxima Beginner’s FAQ

This document is intended to address questions that are asked repeatedly on the Maxima computer algebra system mailing list.

The official Maxima FAQ is located at http://maxima.sourceforge.net/faq.html.

Main contributors are Robert Dodier, Richard Fateman, Stavros Macrakis, Raymond Toy, Barton Willis, Jaime Villate, Andrej Vodopivec and other mailing list subscribers.

Editor: Alexey Beshenov <lyosha@cadadr.org>.

Last modified: Sat, 13 Sep 2008 01:33:57 +0400.

Table of Contents

Algebra

Q: It seems that Maxima could not handle many simple cases of sums of binomial coefficients, infinite sums, etc.

A: Maxima has several different ways of simplifying/evaluating sums:

simplify_sum uses (almost) all algorithms implemented in Maxima to compute the closed form of a sum (simpsum=true, nusum, Zeilberger and also does sums of rational functions using the psi functions). If simplify_sum fails then Maxima probably can’t compute the closed form, so if the user knows about simplify_sum then he does not need to know other functions.

Q: Does Maxima know a way of combining sums which differ only by their bounds using additivity on bounds?

For instance, I don’t know how to simplify (for n>0)

sum(log(i),i,1,n+1)-sum(log(i),i,1,n);

into

log(n+1)

A:

(%i1) sum(log(i),i,1,n+1)-sum(log(i),i,1,n);
                          n + 1           n
                          ====           ====
                          \              \
(%o1)                      >    log(i) -  >    log(i)
                          /              /
                          ====           ====
                          i = 1          i = 1
(%i2) sumcontract(intosum(%));
(%o2)                             log(n + 1)

Q: Could someone explain the concept of the “main polynomial variable” for the Maxima function divide. I do not understand it’s purpose.

A: This is not a Maxima issue, it is a mathematical issue of division with remainder in a polynomial ring, whose coefficients are in a field.

(%i1) divide(x+y, x, x);
(%o1) [1,y]
(%i2) divide(x+y, x, y);
(%o2) [(y+x)/x,0]

Q: How to factor (…) a polynomial over the Gaussian integers?

A: Use gfactor (…).

Q: If my calculations are correct, the roots of

                  6      5      4      3
                 x  + 3 x  + 6 x  + 3 x  + 9 x + 9

should all be expressible by radicals. In fact, they are all polynomials in terms of 2^(1/3) and a cube root of unity.

How can I make Maxima tell me what the roots actually are?

A:

(%i1) p: x^6+3*x^5+6*x^4+3*x^3+9*x+9;
                        6      5      4      3
(%o1)                  x  + 3 x  + 6 x  + 3 x  + 9 x + 9
(%i2) p2: factor(p,q^3-2)$

(%i3) solve(subst(q=2^(1/3),p2),x);
               1/3                    1/3
             (2    + 1) sqrt(3) %i + 2    + 1
(%o3) [x = - --------------------------------,
                            2
      1/3                    1/3              1/3                    1/3
    (2    + 1) sqrt(3) %i - 2    - 1        (2    - 1) sqrt(3) %i + 2    + 1
x = --------------------------------, x = - --------------------------------,
                   2                                       2
      1/3                    1/3                            1/3
    (2    - 1) sqrt(3) %i - 2    - 1        sqrt(3) %i - 2 2    + 1
x = --------------------------------, x = - -----------------------,
                   2                                   2
                    1/3
    sqrt(3) %i + 2 2    - 1
x = -----------------------]
               2

Q: Maxima currently simplifies (x^2)^(1/3) to x^(2/3) which is not correct for all values of x.

(%i1) (x^2)^(1/3);
                                      2/3
(%o1)                                x
(%i2) at (x^(2/3), [x=-1]);
(%o2)                                 - 1

A: Try setting domain to complex:

(%i1) (x^2)^(1/3);
                                      2/3
(%o1)                                x
(%i2) at (%, [x=%i]);
(%o2)                                 - 1
(%i3) domain : complex$

(%i4) (x^2)^(1/3);
                                      2 1/3
(%o4)                               (x )
(%i5) ((-1)^2)^(1/3);
(%o5)                                  1
(%i6) (-1)^(2/3);
                                        2/3
(%o6)                              (- 1)
(%i7) rectform(%);
                                sqrt(3) %i   1
(%o7)                           ---------- - -
                                    2        2

For radcan, sqrt(x^2) is simplified to x, deliberately. radcan does not look at the domain setting. If you don’t want this kind of simplification, don’t use radcan.

Q: radcan(sqrt(x^2-2*x+1)); gives x-1 and not abs(x-1), as expected. Is it a bug in Maxima?

A: it is perhaps inconsistent with this:

(%i1) z : (x-1)^2$
(%i2) sqrt(z);
(%o2)                            abs(x - 1)

though sqrt(expand(z)); does not do that simplification.

Algebraically speaking, sqrt((x-1)^2) has two algebraic values, x-1 and 1-x. note that neither one of them is abs(x-1). radcan chooses one of the values based on which goes to +inf as x goes to +inf.

And if you insist that abs(x-1) is what you want, what consistent theory can you offer to choose a single value for ((x-1)^3)^(1/3)? There are 3 values, 2 which involve 1/2*(sqrt(3)*%i+1), whose powers cover all the roots (that is, they are primitive roots). But x-1 is not a primitive root.

Q: If I make Maxima solve the following example equation it returns results which contain complex parts (%i) but I think the results should be “normal” real numbers.

Example:

(%i1) solve(x^3-4*x+2=0, x);
              sqrt(3) %i   1    - 3/2                 1/3
(%o1) [x = (- ---------- - -) (3      sqrt(37) %i - 1)
                  2        2
           sqrt(3) %i   1
        4 (---------- - -)
               2        2              sqrt(3) %i   1
 + -----------------------------, x = (---------- - -)
       - 3/2                 1/3           2        2
   3 (3      sqrt(37) %i - 1)
                                        sqrt(3) %i   1
                                   4 (- ---------- - -)
   - 3/2                 1/3                2        2
 (3      sqrt(37) %i - 1)    + -----------------------------,
                                   - 3/2                 1/3
                               3 (3      sqrt(37) %i - 1)
      - 3/2                 1/3                 4
x = (3      sqrt(37) %i - 1)    + -----------------------------]
                                      - 3/2                 1/3
                                  3 (3      sqrt(37) %i - 1)

A: How about this:

(%i1) solve(x^3-4*x+2,x)$
(%i2) rectform(%)$
(%i3) trigsimp(%);


(%o3) [x =
                               sqrt(37)                       sqrt(37)
                          atan(---------) - %pi          atan(---------) - %pi
                               3 sqrt(3)                      3 sqrt(3)
          - 2 sqrt(3) sin(---------------------) - 2 cos(---------------------)
                                    3                              3
          ---------------------------------------------------------------------,
                                         sqrt(3)
                       sqrt(37)                       sqrt(37)
                  atan(---------) - %pi          atan(---------) - %pi
                       3 sqrt(3)                      3 sqrt(3)
    2 sqrt(3) sin(---------------------) - 2 cos(---------------------)
                            3                              3
x = -------------------------------------------------------------------,
                                  sqrt(3)
                              sqrt(37)
                         atan(---------) - %pi
                              3 sqrt(3)
                   4 cos(---------------------)
                                   3
               x = ----------------------------]
                             sqrt(3)

All the numbers are real now, but there are trig functions. I think that’s unavoidable.

(%i4) map ('float,%);
(%o4) [x = 0.53918887281089, x = - 2.214319743377535, x = 1.675130870566646]

If you don’t need the exact solution, you can use allroots:

(%i5) allroots(x^3-4*x+2,x);
(%o5) [x = 0.53918887281089, x = 1.675130870566646, x = - 2.214319743377535]

Q: Is there an easy way to get Maxima to say that (1-%i)/(1+%i) = -%i?

A:

(%i1) rectform((1-%i)/(1+%i));
(%o1)                               - %i
(%i2) rat((1-%i)/(1+%i)), algebraic;
(%o2)/R/                            - %i

Q: In Maxima, the polynomial (x-4)^2*(x-2)^3*(x-1)*(x+3) has two complex roots and only one integer root. Try it yourself:

(%i1) (x-4)^2*(x-2)^3*(x-1)*(x+3)$
(%i2) allroots(%);
(%o12) [x = 1.0, x = 4.213716426528219E-5 %i + 1.999957682213045,
x = 1.999957682213045 - 4.213716426528219E-5 %i, x = 2.000084638451563,
x = 3.999949310178406, x = - 2.999999999941271, x = 4.000050686885212

A: allroots is a numerical solver and roots are approximate. See that imaginary parts are near zero.

An alternative might be to use polyroots, which you can get from load(jtroot3). This is essentially the same algorithm as allroots, but supports bfloat arithmetic. If I set fpprec to 32, I get the roots:

[9.9999999999999999999999999999999b-1 - 4.2370458776519188767231514736868b-33 %i,
2.0967098044841547417705251129837b-11 %i + 1.999999999996426866873956584756b0,
1.9509670336642699191457056558385b-31 %i - 3.0000000000000000000000000000005b0,
1.9999999999836285277997026851267b0 - 1.3577973543378557813146327177067b-11 %i,
2.0000000000199446053263407300922b0 - 7.3891245014629895950001484527873b-12 %i,
4.8135821109982436660981576396279b-15 %i + 4.0000000000000008869298470996269b0,
3.9999999999999991130701529003988b0 - 4.8135821109982534157333151111319b-15 %i]

And prt_root_error tells me that the estimated errors in the roots are:

[8.234297913969850316964397927954b-31,
4.9446128309052513207164850565785b-9,
6.0819981869720905951308760004527b-31,
4.9736209677990489185125295584204b-9,
4.8869211713208376093083688645524b-9,
3.5862126603415141455614873576955b-15,
3.57319044991954585353389212079b-15]

Thus, the single root at 1 is very accurate, as is the single root at −3. The rest have reduced accuracy because they are repeated.

Q: is there a way to solve a complex equation with conjugate as:

(%i1) solve([conjugate(z)-3*%i*z-3+6*%i=0], [z]);
                                     6 %i - 3
(%o1)                           [z = --------]
                                     3 %i - 1

wich gives a wrong solution

                                     6 %i - 3
(%o1)                           [z = --------]
                                     3 %i - 1

A: Maxima assumes that variables are real-valued by default (even with domain:complex), so conjugate(z) will yield z. This can be solved by declare(z,complex).

But alas, Maxima’s solve command knows nothing about conjugate, so this doesn’t help with the original problem.

The easiest way to solve it is to explicitly solve for the real and imaginary parts:

(%i1) z:a+%i*b$
(%i2) conjugate(z)-3*%i*z-3+6*%i;
(%o2)               - 3 %i (%i b + a) - %i b + a + 6 %i - 3
(%i3) [realpart(%),imagpart(%)];
(%o3)                    [3 b + a - 3, - b - 3 a + 6]
(%i4) solve(%,[a,b]);
                                     15      3
(%o4)                          [[a = --, b = -]]
                                     8       8
(%i5) subst(%[1],z);
                                   3 %i   15
(%o5)                              ---- + --
                                    8     8

Q: algsys([2*z=1, z*z+c=z], [c]); does not work.

A:

algsys([2*z=1, z*z+c=z], [c,z]);

Q: I don’t understand the following results from solve.

(%i1) solve([y = 2*x, y = sqrt(x^2+3)],[x,y]);
algsys' cannot solve - system too complicated.
 -- an error.  To debug this try debugmode(true);

A: In general solve cannot solve expressions which contain sqrt(foo(x)). to_poly_solve helps:

(%i1) load(topoly_solver)$
(%i2) to_poly_solve(2*x=sqrt(x^2+3),x);
(%o2)                              [[x = 1]]

Q: Is there a recommended way to do the composition of functions?

A: If you deal with expressions:

(%i1) comp(g,f,x) := at(g,x=f)$
(%i2) comp(sqrt(1-x^2),sqrt(x),x);
(%o2)                             sqrt(1 - x)
(%i3) comp(asin(x),sin(x),x);
(%o3)                                  x

If you deal with functions:

(%i1) comp(g,f,x) := g(f(x))$
(%i2) comp(lambda([x],sqrt(1-x^2)),'sqrt,x);
(%o2)                             sqrt(1 - x)
(%i3) comp('asin,'sin,x);
(%o3)                                  x
(%i4) comp(g,f,[x]) := g(apply(f,x))$
(%i5) comp(lambda([x],x^2),'mod,24,5);
(%o5)                                 16

Matrices and Linear Algebra

Q: Is there even faster way to solve linear equations than solve or matrix inversion?

A: See fast_linsolve.

Q: I have a matrix A (for example 3×3), and a vector X=matrix([x],[y],[z]). What are the instructions for solving the linear system A.X=0?

A: For a invertible coefficient matrix, you can use the (undocumented) function linsolve_by_lu; for other cases, I think you’ll have to convert to equations and use linsolve (or solve or algsys).

Examples:

Solve using rational numbers—doesn’t compute the matrix condition number:

(%i1) linsolve_by_lu(matrix([5,7],[9,6]), matrix([1],[1]));

                                 [ 1  ]
                                 [ -- ]
                                 [ 33 ]
(%o1)                           [[    ], false]
                                 [ 4  ]
                                 [ -- ]
                                 [ 33 ]

Solve using IEEE doubles—returns an upper bound for the matrix condition number as the second list entry:

(%i2) linsolve_by_lu(matrix([5,7],[9,6]), matrix([1],[1]), 'floatfield);

                    [ 0.03030303030303 ]
(%o2)              [[                  ], 6.835016835016835]
                    [ 0.12121212121212 ]

Multiple right-hand-sides are okay—find the inverse:

(%i3) linsolve_by_lu(matrix([5,7],[9,6]), matrix([1,0],[0,1]));
                             [   2    7   ]
                             [ - --   --  ]
                             [   11   33  ]
(%o3)                       [[            ], false]
                             [  3      5  ]
                             [  --   - -- ]
                             [  11     33 ]

Solve using bigfloats:

(%i4) linsolve_by_lu(matrix([5,7],[9,6]), matrix([1],[1]), 'bigfloatfield), fpprec : 20;

                  [ 3.03030303030303b-2  ]
(%o4)            [[                      ], 6.835016835016835]
                  [ 1.212121212121212b-1 ]

Q: I am using linsolve() for the solution of linear equations:

(%i1) ysol : linsolve([a11*x1+a12*x2+a13*x3],[x1,x2,x3]);
                      %r1 a13 + %r2 a12
(%o1)         [x1 = - -----------------, x2 = %r2, x3 = %r1]
                             a11

I want to transform the obtained solution to the form xsol = A . r.

A: Try this:

(%i1) [x1 = -(%r1*a13+%r2*a12)/a11,x2 = %r2,x3 = %r1]$

(%i2) map('rhs,%);
                       - %r1 a13 - %r2 a12
(%o2)                 [-------------------, %r2, %r1]
                               a11
(%i3) coefmatrix(%,[%r1,%r2]);
                             [   a13    a12 ]
                             [ - ---  - --- ]
                             [   a11    a11 ]
(%o3)                        [              ]
                             [   0      1   ]
                             [              ]
                             [   1      0   ]

Why map rhs on to the equations? Consider:

(%i4) coefmatrix([x+y=z, x-y=q],[x,y]);
                                [ 1   1  ]
(%o4)                           [        ]
                                [ 1  - 1 ]
(%i5) augcoefmatrix([x+y=z, x-y=q],[x,y]);
                              [ 1   1   - z ]
(%o5)                         [             ]
                              [ 1  - 1  - q ]

Q: How to apply the modulo function to matrices?

A: mod is the modulo function. If you want to apply mod 26 to a matrix m, use matrixmap:

(%i1) m: matrix ([15,2,7], [8,10,23], [0,2,8])$

(%i2) matrixmap(lambda([x],mod(x,26)),m);
                               [ 15  2   7  ]
                               [            ]
(%o2)                          [ 8   10  23 ]
                               [            ]
                               [ 0   2   8  ]

If you don’t like the lambda construction you can always create a little function first:

(%i1) mod26(x):=mod(x,26)$
(%i2) matrixmap(mod26,m);

Another way to perform modular arithmetic is to set a global modulus, which applies to all calculations on Canonical Rational Expressions, which you create using rat(…). For example:

(%i1) m: matrix([15,2,7],[8,10,23],[0,2,8])$

(%i2) determinant(m);    /* calculate over Z */
(%o2)                                 494
(%i3) modulus:23;
(%o3)                                 23
(%i4) rm:rat(m);
                                [ - 8  2   7 ]
                                [            ]
(%o4)/R/                        [  8   10  0 ]
                                [            ]
                                [  0   2   8 ]
(%i5) determinant(rm);    /* calculate over Z23 */
(%o5)/R/                              11
(%i6) rm^^-1;    /* inverse over Z23 */
                               [  1     4    2 ]
                               [               ]
(%o6)                          [ - 10  - 10  3 ]
                               [               ]
                               [ - 9   - 9   8 ]
(%i7) rm . rm^^-1;
                                  [ 1  0  0 ]
                                  [         ]
(%o7)/R/                          [ 0  1  0 ]
                                  [         ]
                                  [ 0  0  1 ]

Q:

(%i1) func(n,m):=for i:1 thru n do for j:1 thru m do A[i,j]:i+j$
(%i2) A;
(%o2)                                 A
(%i3) A[1,2];
(%o3)                                 3

But what is A? Matrix? Array?

I tried to define matrix with array’s help:

(%i4) func(n,m):=(array(A,flonum,n,m),for i:1 thru n do for j:1 thru m
          do A[i,j]:i+j,M:genmatrix(A,n,m))$
(%i5) func(4,4);

But without success :-(

Element and array type do not match:
1234
#0: func(n=4,m=4)
 -- an error.  To debug this try debugmode(true);

I need to define all elements of matrix in function. How can I do this?

A: You have found one of the more confusing areas of Maxima, I’m afraid.

The three different tests you did in fact operate on three different types of object in Maxima: sparse arrays (assigning a[..] with no declaration), matrices (matrix/genmatrix/etc.), and dense arrays (explicit array declaration).

First of all, I would recommend that you not use dense arrays at all. Those are best reserved for optimization of performance by advanced users.

Sparse arrays allow you to store values very flexibly, but do not work as single assignable objects. For example:

(%i1) qq[3]: 5$

(%i2) qq["my dog"]: 'fido $

(%i3) qq[-2/3] : "minus two thirds";$

(%i4) qq[3];
(%o4)                                  5

But you can’t do much with qq by itself—you can’t print it, you can’t perform arithmetic on it (e.g. add elements index-by-index with another sparse array), etc. (but see arrayinfo and listarray for more advanced operations).

You can write matrices explicitly, e.g. matrix([1,0],[0,1]).

You can create matrices using genmatrix:

(%i5) genmatrix( lambda([i,j], i+j), 4,4);
                                [ 2  3  4  5 ]
                                [            ]
                                [ 3  4  5  6 ]
(%o5)                           [            ]
                                [ 4  5  6  7 ]
                                [            ]
                                [ 5  6  7  8 ]

To make a matrix from a sparse array, use genmatrix:

(%i6) for i thru 4 do for j thru 4 do A[i,j]: i+j$

(%i7) genmatrix(A,4,4);
                                [ 2  3  4  5 ]
                                [            ]
                                [ 3  4  5  6 ]
(%o7)                           [            ]
                                [ 4  5  6  7 ]
                                [            ]
                                [ 5  6  7  8 ]

You can also modify the elements of a matrix:

(%i8) MA: genmatrix(A,3,3)$
(%i9) for i thru 4 do MA[i,i]: 0$
(%i10) MA;
                                 [ 0  3  4  5 ]
                                 [            ]
                                 [ 3  0  5  6 ]
(%o10)                           [            ]
                                 [ 4  5  0  7 ]
                                 [            ]
                                 [ 5  6  7  0 ]

There are various other constructors for matrices, e.g.

(%i11) ident(2);
                                   [ 1  0 ]
(%o11)                             [      ]
                                   [ 0  1 ]

Try working with these constructs.

Q: How can I build a matrix with all its elements zero with a simple command?

A: You are looking for zeromatrix:

(%i1) zeromatrix(4,4);

More generally, genmatrix will construct a matrix using any function of the indices you like:

(%i2) genmatrix(lambda([i,j],0),4,4);

Q: Is the following behavior desired?

(%i1) M: matrix([1,2,3])$
(%i2) M . transpose(M);
(%o2)                                 14

One would expect multiplying a 1×n matrix by a n×1 matrix to return a 1×1 matrix rather than a scalar.

A:

(%i1) scalarmatrixp : false$
(%i2) M : matrix([1,2,3])$
(%i3) M . transpose(M);
(%o3)                               [ 14 ]

Q: charpoly very, very slow

(%i1) matrix([0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4],
             [1,0,2,1,2,1,3,2,2,1,3,2,3,2,4,3],
             [1,2,0,1,2,3,1,2,2,3,1,2,3,4,2,3],
             [2,1,1,0,3,2,2,1,3,2,2,1,4,3,3,2],
             [1,2,2,3,0,1,1,2,2,3,3,4,1,2,2,3],
             [2,1,3,2,1,0,2,1,3,2,4,3,2,1,3,2],
             [2,3,1,2,1,2,0,1,3,4,2,3,2,3,1,2],
             [3,2,2,1,2,1,1,0,4,3,3,2,3,2,2,1],
             [1,2,2,3,2,3,3,4,0,1,1,2,1,2,2,3],
             [2,1,3,2,3,2,4,3,1,0,2,1,2,1,3,2],
             [2,3,1,2,3,4,2,3,1,2,0,1,2,3,1,2],
             [3,2,2,1,4,3,3,2,2,1,1,0,3,2,2,1],
             [2,3,3,4,1,2,2,3,1,2,2,3,0,1,1,2],
             [3,2,4,3,2,1,3,2,2,1,3,2,1,0,2,1],
             [3,4,2,3,2,3,1,2,2,3,1,2,1,2,0,1],
             [4,3,3,2,3,2,2,1,3,2,2,1,2,1,1,0])
(%i2) charpoly(%,x);

A: charpoly computes the determinant using a slow algorithm. If you set ratmx to true, a faster algorithm is used.

Analysis and Statistics

Q: I want to teach Maxima how to properly treat the derivative of the conjugate function.

With the end result that I will get the behavior:

(%i1) diff (conjugate (f(x)), x);
                                      d
(%o1)                       conjugate(-- (f(x)))
                                      dx

A: Here is such a rule.

(%i1) declare (e, complex)$
(%i2) matchdeclare (e, all, u, mapatom)$
(%i3) tellsimpafter ('diff (conjugate (e), u), conjugate (diff (e, u)))$

This depends rather sensitively on Maxima’s assumptions about how diff and conjugate work; e.g. there must be a quote mark on diff but not on conjugate. I hope we can eventually regularize some things so that making up basic rules like this doesn’t depend on peculiar factoids… Anyway here are a couple of examples of the rule in action.

(%i4) diff (conjugate (f(x)), x);
                                      d
(%o4)                       conjugate(-- (f(x)))
                                      dx
(%i5) diff (conjugate (f(x) + g(x)), x);
                           d                      d
(%o5)            conjugate(-- (g(x))) + conjugate(-- (f(x)))
                           dx                     dx

Q: I’m having trouble getting a Taylor series expansion of an expression involving the derivative.

fo:taylor('diff(u(t),t)+a*u(t),t,0,4);

Gives me the expansion of f(t) around t=0 and I can work with that. Then I am stumped with how to get the Taylor series expansion of f(2*t) around t=0. Something like

fo:taylor('diff(u(2*t),t)+a*u(2*t),t,0,2)

really doesn’t work.

A: Maybe pdiff will do what you want:

(%i1) load(pdiff)$

(%i2) diff(u(t),t)+a*u(t);
(%o2)                          a u(t) + u (t)
                                         1
(%i3) subst(t = t+2*dt,%);
(%o3)                   a u(t + 2 dt) + u (t + 2 dt)
                                         1

(%i4) taylor(%,dt,0,2);
                 2
(%o4) a u(t) + dt  (2 u (t) + 2 a u (t)) + dt (2 u (t) + 2 a u (t)) + u (t) + ...
                       3           2              2           1        1

Q: Does Maxima know how to express special integrals in terms of special functions? I suggest developers to modify integrate() to use Ei(t), Si(t), Ci(t) in output.

A: Adding new cases to the integration routines is not all that easy.

To do this right you would need to determine if an answer can be expressed as some function f(Ei(g(x)), Si(h(x)), …)

For example, would you know how to recognize that an integral is Ei(sqrt(x^2+log(x))? If the check is only going to work for exactly the derivative of Ei(x) it is not very useful.

In fact, if that is all you want, you can use

tellsimp(integrate(exp(x)/x,x), …)

and you don’t need any developers at all.

Q: When I try to integrate some functions (part of a double integral) and I fill in some bounds, Maxima always keeps asking: is foo positive, negative or zero.

Can this question turned off?

A: In general, no. However if Maxima is asking about the sign of a symbol (as opposed to a more general expression) then you can set a default response via assume_pos and assume_pos_pred.

Q: I want to evaluate the improper integral

(%i1) integrate(1/(a+sin(x)),x,0,2*%pi);
Is  (a - 1) (a + 1)  positive, negative, or zero?

p;
          2
Is  sqrt(a  - 1) + a  positive or negative?

p;
    !      2         !
Is  !sqrt(a  - 1) - a! - 1  positive, negative, or zero?

n;    /* I had to do a slight calculation */
          2
Is  sqrt(a  - 1) + a - 1  positive, negative, or zero?

p;
                                      2
                          2 %pi sqrt(a  - 1) - 2 %pi a
(%o1)                   - ----------------------------
                                    2         2
                            a sqrt(a  - 1) - a  + 1

The answer is correct, but I wonder what kind of procedure was followed and I had to specify all those values.

I tried to replace a with b+1 with b positive, but I received 4 questions again with only the first actually needed.

The issue with questions of that type is something that is planned to be fixed? Where can I read details about this, so I don’t post trivial questions myself?

A: You could start with assume(a>1). This gets rid of the first two questions.

The first two questions come about because Maxima is trying to see if the integrand has any poles on the real axis on the interval from 0 to 2*%pi. (I discovered this by pressing C-c and looking at the backtrace.)

After that, Maxima is applying the substitution exp(%i*x) = y and doing a contour integral around the unit circle. This is done via residues, so the questions come from computing the poles of the new integrand.

Definite integration is pretty much contained in the Lisp file src/defint.lisp. There are some comments in the file describing what’s happening.

Q: Sometimes Maxima asks questions about your variables whether they are greater than zero or etc. Then, it answers your question as your data. Is there a way to make Maxima give all the results for different values and combinations of variables at the same time? If yes how? If no, would not it be a good feature to implement?

A: It’s a great idea, trhere is an add-on package share/contrib/noninteractive which intercepts asksign questions and turns them into conditional expressions. However noninteractive sometimes fails badly.

load (noninteractive);
integrate (x^k, x);
 => if equal(k+1,0) then log(x) else x^(k+1)/(k+1)

See the comments in noninteractive.mac for other examples.

Q: How do I find the Laplace transform and inverse Laplace transform of some function?

A: See laplace(f(t),t,s) and ilt(F(s),s,t).

(%i1) laplace (exp(t)*cos(t),t,s);
                                    s - 1
(%o1)                            ------------
                                  2
                                 s  - 2 s + 2
(%i2) ilt(%,s,t);
                                    t
(%o2)                             %e  cos(t)
(%i3) laplace ('integrate(f(x)*g(t-x), x, 0, t), t, s);
(%o3)               laplace(f(t), t, s) laplace(g(t), t, s)
(%i4) laplace (diff(f(t),t,3), t, s);
                  !
         2        !                    !
        d         !           d        !          3
(%o4) - --- (f(t))!      - s (-- (f(t))!     ) + s  laplace(f(t), t, s)
          2       !           dt       !
        dt        !                    !t = 0
                  !t = 0
                                                                              2
                                                                      - f(0) s

Q: desolve cannot solve some differential equations that ode2 can.

A: desolve finds the Laplace transform of each equation and, only in the simple cases when the resulting system is an algebraic system, it solves it.

ode2 solves various different types of first and second order differential equations.

Q:

(%i1) h : -(k^(k/(k-1))*(2*k-1)*2^(k/(k-1)))/
               ((k-1)*(k*2^((2*k)/(k-1))-2^((2*k)/(k-1))-2*k+1))$

(%i2) limit(h,k,1/3,plus);
                                      log(3)
                                      ------
                                        2
                                  3 %e
(%o2)                             -----------
                                  160 sqrt(2)
(%i3) limit(h,k,1/3,minus);
                                      log(3)
                                      ------
                                        2
                                  3 %e
(%o3)                             -----------
                                  160 sqrt(2)

Both limits are wrong: the first one is +∞ and the second one is −∞. Can one really trust Maxima regarding the calculation of limits?

A: It has been suggested that the limit calculations should be reprogrammed along the lines of a recent PhD dissertation by Dominik Gruntz.

Your example answers your own question about whether you can really trust the current program :-)

Q: Could Maxima compute the limit

(%i1) limit (gamma ((x+1)/2)/(sqrt(x)*gamma (x/2)),x,inf);
                                        x   1
                                       (- - -)!
                                        2   2
(%o1)                     limit    ----------------
                          x -> inf  x
                                   (- - 1)! sqrt(x)
                                    2

(I believe that the answer should be 1/sqrt(2))?

A:

(%i1) gamma ((x+1)/2)/(sqrt(x)*gamma (x/2))$
(%i2) tlimit(%, x, inf);
                                       1
(%o2)                               -------
                                    sqrt(2)

Q: Can I use Maxima to fit some data set to my desired function?

A: There is an add-on package lsquares by which more general functions can be fit to data. First one constructs an expression for the mean-square error, then one finds parameter estimates by minimizing that error.

?? lsquares

finds the documentation for lsquares, with several examples.

Q: Are there some way to optimize a function (non-linear) constrained to some equations and/or inequalities?

A: There is augmented_lagrangian_method which can solve problems with equality constraints on a nonlinear objective function. For inequality constraints, there isn’t anything in Maxima to directly solve such problems. Maybe you can try to rephrase the problem as an unconstrained problem, e.g. a constraint x > 0 could be translated into an unconstrained variable y = log(x).

Or maybe the objective function is close enough to linear so you could apply the simplex method, which Maxima has.

augmented_lagrangian_method is not the strongest method for problems with equality constraints.

Q:

(%i1) solve(sin(x) + x = 2.5*x + 1.5, x);

                                  2 sin(x) - 3
(%o1)                        [x = ------------]
                                        3

Isn’t it possible to solve this equation numerically?

A: You can plot the stuff returned by solve, and use that to guess a range or ranges for find_root. (find_root only tries to find one root.)

(%i1) plot2d ([(2*sin(x)-3)/3, x], [x, -%pi, %pi]);
(%i1) find_root(sin(x)+x-2.5*x - 1.5, x, -%pi, %pi);
(%o1)                        - 1.663786360949127

Q: I’m trying to find all points where two functions intersect, using Maxima:

(%i1) log2(x) := log(x)/log(2)$
(%i2) eq1 : y = 64*n*log2(n)$
(%i3) eq2 : y = 8*n^2$
(%i4) solve([eq2,eq1],[n,y]);

`algsys' cannot solve - system too complicated.
 -- an error.  Quitting.  To debug this try debugmode(true);

A: Maxima cannot solve this system of equations symbolically because it is not solvable in closed form. Using the eliminate function, you can reduce it to an equation in n only:

(%i1) eq12: eliminate([eq1,eq2],[y]);
(%o1)                     [8 n (8 log(n) - log(2) n)]
(%i2) sols: solve(eq12,n);
                                         8 log(n)
(%o2)                        [n = 0, n = --------]
                                          log(2)

The first solution is correct (despite the log(n) term…); the second is a transcendental equation and can’t be solved symbolically using the usual elementary functions. At this point, you can use find_root:

(%i3) find_root(sols[2],n,.1,10);
(%o3)                                1.1
(%i4) find_root(sols[2],n,10,100);
(%o4)                               43.56

Q: Help me understand this inconsistency?

(%i1) limit(2^x,x,inf,plus);
(%o1) inf
(%i2) limit(2^x,x,inf,minus);
(%o2) inf – why isn’t this ‘-inf’ ?

A: The fourth argument is only used for finite limits.

Q: How to set Maxima in degree mode for trigonometric functions?

A: Maxima trigonometric routines take arguments expressed in radians. There is no way to tell Maxima to use degrees instead. However, you can easily define trigonometric routines with degree arguments:

(%i1) deg_to_rad(d) := %pi*d/180$
(%i2) rad_to_deg(r):= r*180/%pi$
(%i3) sind(x):=sin(deg_to_rad(x))$
(%i4) asind(x) := rad_to_deg(asin(x))$
...
(%i5) sind(45);
                                       1
(%o5)                               -------
                                    sqrt(2)
(%i6) sind(30);
                                       1
(%o6)                                  -
                                       2
(%i7) sind(15);
                                      %pi
(%o7)                             sin(---)
                                      12
(%i8) ev (%, numer);
(%o8)                         .2588190451025207

Q: Is there a way to do regressions and correlations with Maxima?

A: Take a look at the lsquares and stats packages.

You might consider working with R which is a general statistical package.

Q: How can I solve an inequality in Maxima?

A: You can try fourier_elim. There is a pre-processor to fourier_elim that allows it to solve some nonlinear inequalities.

(%i1) load(fourier_elim)$

(%i2) fourier_elim([abs(x - abs(5-x)) < 1],[x]);
(%o2)                          [2<x,x<3]

(%i3) fourier_elim([x + y < 1, x - y > 9],[x,y]);
(%o3)                     [y+9<x,x<1-y,y<-4]

(%i4) fourier_elim([max(-x,x) < 7 * x + 1],[x]);
(%o4)                           [-1/8<x]

Plotting

Q: I cannot obtain the 3D plotting of two functions (and more functions too):

plot3d([7*x-3*y, x^2-y^2], [x,-5,5], [y,-5,5])$

A: plot3d can’t plot multiple surfaces. You can use the draw package.

There are some examples here: http://www.telefonica.net/web2/biomates/maxima/gpdraw

Q: I have a function defined as follows:

F(x)=(e^x)/3 for x<=0
F(x)=1 for x>0

I need to plot this function. How do I do that bearing in mind that the function definition changes according to the values of x?

A:

You can use an if-then-else construct to define your function.

(%i1) f(x):= if x<=0 then (%e^x)/3 else 1$
(%i2) plot2d(f(x),[x,-2,2],[y,0,1.1])$

The plot shows a vertical connection between the parts. If you don’t want to see this, define the two parts separately. In each I leave one part undefined, which remains unplotted.

(%i3) f1(x):= if x<=0 then (%e^x)/3$
(%i4) f2(x):= if x>0 then 1$
(%i5) plot2d([f1(x),f2(x)], [x,-2,2], [y,0,1.1])$

With the draw package:

load(draw)$
draw2d(explicit (%e^x/3, x, -3, 0),
       explicit (1,      x,  0, 3),
       yrange     = [-2,2],
       grid       = true,
       point_type = filled_circle,
       point_size = 2,
       points ([[0,1/3]]))$

Q: How do I get color in an eps file using draw2d’s image function?

A: Use terminal = eps_color

Notations, Evaluation, Simplification

Q: Numeric expressions are evaluated, but that doesn’t work for boolean expression:

(%i1) 5>3;
(%o1)                                5 > 3

A:

(%i1) is(5>3);
(%o1)                               true

The is command is not necessary in other contexts:

(%i2) if 5>3 then 1 else 0;
(%o2)                                  1
(%i3) 5>3 or 2>4;
(%o3)                                true

Q: I have two expressions, test1=1/sqrt(2) and test2=2*sqrt(2-sqrt(2))/sqrt(2), which I would like Maxima to display as sqrt(2)/2 and sqrt(2*(sqrt(2)-1)), respectively. Both involve rationalizing the denominator, and the second involves pulling terms inside a radical (a la rootscontract) and factoring inside. Is there a straightforward way to obtain these results?

A: Try ratsimp with algebraic : true

(%i1) algebraic : true$
(%i2) ratsimp(test1=1/sqrt(2));
(%o2) test1=sqrt(2)/2
(%i3) ratsimp(test2=2*sqrt(2-sqrt(2))/sqrt(2));
(%o3) test2=sqrt(2-sqrt(2))*sqrt(2)

Q: I want Maxima to convert exponential expression to hyperbolic one. How can I do that? I tried demoivre but it converts exponential to trigonometric but not to hyperbolic.

A: Here’s the trick: change the exponents to imaginary, use demoivre, then change back. For example:

(%i1) exp(3*x)/(exp(x)-exp(2*x));
                                       3 x
                                     %e
(%o1)                             -----------
                                    x     2 x
                                  %e  - %e
(%i2) subst(%i*x,x,%);
                                     3 %i x
                                   %e
(%o2)                          -----------------
                                 %i x     2 %i x
                               %e     - %e
(%i3) demoivre(%);
                            %i sin(3 x) + cos(3 x)
(%o3)            ---------------------------------------------
                 - %i sin(2 x) - cos(2 x) + %i sin(x) + cos(x)
(%i4) subst(x/%i,x,%);
                             sinh(3 x) + cosh(3 x)
(%o4)             -------------------------------------------
                  - sinh(2 x) - cosh(2 x) + sinh(x) + cosh(x)

Q: Why do Maxima change sequence of variables?

(%i1) (a + b)^2;
                                          2
(%o1)                              (b + a)

A: There are two things going on here. One is that Maxima tries to arrange expressions into a so-called canonical form. The other is that Maxima has ideas about how to display the canonical form of an expression.

You can change this behavior via the global flag powerdisp.

(%i1) a + b + c;
(%o1)                              c + b + a
(%i2) powerdisp : true$

(%i3) a + b + c;
(%o3)                             a + b + c

To judge by the name, powerdisp was invented to control the display of polynomials. But it actually affects all + expressions.

Q: As a Maxima newbie I’ve trouble understanding how to achieve a certain transformation.

Especially I would like the following term:

  2 vx z - r z - l z - 2 n (x - vx) - 2 vz x + r vz + l vz
- --------------------------------------------------------
                           r - l

to be written as:

2(n + vz)      -2 vx    r + l       2 n vx + r vz + l vz
--------- x + (----- +  -----) z +  --------------------
  r - l        r - l    r - l              r - l

A: There is actually no reason to expect Maxima to “simplify” an expression to some arbitrary form that you happen to like; Maxima has its own particular repertoire of “kinds of simplifications” and getting exactly the ordering you want may be contrary to its rules.

On the other hand, there is a package that tries to come up with a formatting arrangement that might allow for something like this form. Bruce Miller wrote it, and it is called format. See http://math.nist.gov/~BMiller/computer-algebra/

Q: I have a matrix with very long trigonometric expressions. Now when I want to get the jacobian of that matrix there is an error message that the expression is too long to be shown.

A: Try setting display2d : false and something like

print_matrix(M) := (
    for i : 1 thru length(M) do
        for j : 1 thru length(M[i]) do
            print (concat('M,"[",i,"]","[",j,"] :"), M[i,j])
)$

Exporting to TeX and rendering the TeX expression may be useful too.

Q: Using trigreduce(y) leads to terms of cos(w1), cos(w2), cos(2*w1+w2) … and so on. I’d like to sort these terms according to those cos() coefficients.

A: The sort function sorts all kinds of expressions.

sort ([cos(w1 + w2), cos(w2^2), cos(w1), cos(2*w1 + w2), cos(w2)]);
 => [cos(w1), cos(w2), cos(w2 + w1), cos(w2 + 2 w1), cos(w2^2)]

Maxima has its own idea about how to order non-atomic expressions. If you don’t like that order, you can tell sort to use a specific comparison function.

? sort

says more about that.

Q: I have a function, f(x), that when Taylor expanded around x=1 in Maxima gives:

                           2
         M_tot N   (M_tot N  + 2 M_tot N) (x - 1)
(%o1)/T/ ------- + ------------------------------
            2                    12
                                                2                     2
                                        (M_tot N  + 2 M_tot N) (x - 1)
                                      - -------------------------------
                                                      24
+ . . .

Is there a flag that will simplify the coefficients of (x-1) separately?

A: Try things like map('factor, %) or map('ratsimp, %), or …

Floating Point

Q: When I type something like sqrt(2) or 3/2 Maxima shows it exactly as I wrote it. Is there a way I can make it show the results in decimal form?

A: Yes, use float(…) or ev(…, numer).

Maxima doesn’t do this by default because Maxima keeps everything in exact form until you tell it you want an approximate numeric computation. This permits its results to be more accurate.

Q: Can someone please explain what the distinction is between floats and bigfloats, and why one would select one over the other?

A: bfloat is an arbitrary precision type with fpprec significant digits:

(%i1) float(%pi);
(%o1)                          3.141592653589793
(%i2) bfloat(%pi);
(%o2)                         3.141592653589793b0
(%i3) fpprec : 100;
(%o3)                                 100
(%i4) float(%pi);
(%o4)                          3.141592653589793
(%i5) bfloat(%pi);
(%o5) 3.1415926535897932384626433832795028841971693993751058209749445923078164\
06286208998628034825342117068b0

See the Maxima Manual, section 10 (“Floating Point”).

Q: How to change the calculational precision locally?

A: It’s very straightforward:

block([fpprec:100], bfloat(%pi))

Changing the precision of reading literal constants is a little harder. There are three approaches:

  1. use ''(fpprec: …) to change the precision before parsing the rest of the expression. I find this very ugly and dangerous (like many uses of ''). And though you can do this “locally” in the sense of limiting it to one high-level input, you have no more control than that.
  2. use exact notation instead of floating notation for your numbers, e.g. block([fpprec:100], bfloat(sqrt(1000000000000000000001*10^-21)));
  3. use eval_string or parse_string.

I find (2) the most elegant, but perhaps a little clumsy notationally.

Q: Are there some functions for rounding/truncating?

ceilingnext largest integer
floornext smallest integer
truncatenext integer closer to 0
roundnearest integer, preferring evens to odds for ties

Q: How can I calculate machine epsilon in Maxima?

A: What you want are constants available in any ANSI Common Lisp

(%i1) :lisp least-positive-double-float

4.9406564584124654E-324
(%i1) :lisp  double-float-epsilon

1.1107651257113995E-16
(%i1) :lisp most-positive-double-float

1.7976931348623157E308

This should be the same for any IEEE double float implementation, though it may depend upon whether you allow only normalized numbers.

(%i1) :lisp least-positive-normalized-double-float

2.2250738585072014E-308
(%i1) :lisp least-positive-normalized-single-float

2.2250738585072014E-308

You can compute these items by fairly obvious simple loops, assuming the representation is radix 2, though. I suggest you use (expt 2.0d0 n) rather than (expt 2 n) which is not floating point at all.

Q:

for i: -1 step 0.1 thru 2 do display(i);
                                   i = - 1
                                   i = - 0.9
                                   i = - 0.8
                                   i = - 0.7
                                   i = - 0.6
                                   i = - 0.5
                                   i = - 0.4
                                   i = - 0.3
                                   i = - 0.2
                                   i = - 0.1
                         i = - 1.3877787807814457E-16
                                    i = 0.1
                                    i = 0.2
                                    i = 0.3
                                    i = 0.4
                                    i = 0.5
                                    i = 0.6
                                    i = 1.7
                             i = 1.800000000000001
                             i = 1.900000000000001

Why do steps at 1.8 and 1.9 are displayed with 16 digits, unlike the other steps?

A: Because 0.1 cannot be represented exactly as a floating-point number, and repeatedly adding 0.1 will introduce (more) roundoff errors. This is also why Maxima prints out i = - 1.3877787807814457E-16 instead of 0. Extra digits are needed because they are needed to represent the number.

Maxima has rational numbers. Why not use them? Or count from -1 to 20 and divide i by 10.0. The results will be more accurate, but will still have roundoff issues.

Input and Output

Q: How and where I can save a Maxima session?

A: The function save copies functions and variables in the current session into a file. load restores the stuff saved. At the input prompt,

? save

and

? load

should display some information about those functions.

writefile copies stuff from the console to a file, so all the input and output is stored. The content of the file created by writefile cannot be restored.

Q: As an old Mathcad user I’m used to save some of my most used formulas as a worksheet. Is there some similar functionality in Maxima? Basically, when I’ve defined a function, is it possible to save that function somewhere, close Maxima, open it again, and use that function?

A: Probably the easiest thing is grind(myfunction) which prints the function definition, then copy & paste that into a file. Then next time you want to use that function, load("myfile.mac") loads it. There are other ways to achieve that effect, but that is simplest.

Q: Does maxima have the functionality to import/export data table/matrices from .csv files?

A: Try numericalio.

Q: How to make Maxima output the result as a single-line expression?

A: To convert to a 1-line string (which you can’t compute with):

string(…);

To have Maxima give all its output in linear form:

display2d:false$

Q: When exporting to TeX, cos(x)*y is printed as \cos x\,y and my students could read this as cos(x*y).

A: Try this:

:lisp (defun $detex (sym) (remprop sym 'tex) (remprop sym 'texrpb))
map('detex,[sin,cos,tan,cot])$

This should make everything output as a regular function application.

Q:

(%i1) matrix([1,1,1],[1,1,1],[1,1,1]);
                                  [ 1  1  1 ]
                                  [         ]
(%o1)                             [ 1  1  1 ]
                                  [         ]
                                  [ 1  1  1 ]
(%i2) tex(%)$
$$\pmatrix{1&1&1\cr 1&1&1\cr 1&1&1\cr }$$

\pmatrix does not render properly for MimeTeX, which is the program we’re using to display the output of Maxima.

However, this would work:

$$\begin{pmatrix}1&1&1\\ 1&1&1\\ 1&1&1\end{pmatrix}$$

A: This is LaTeX, not plain TeX. Maxima traditionally generated plain TeX.

load("mactex-utilities") loads an alternate TeX function for matrix which outputs \begin{pmatrix} ... \end{pmatrix}. It also changes the TeX output for quotients. If you want just the matrix stuff, I guess you can copy & paste just that part to maxima-init.lisp or something.

Also it may be easier to define \pmatrix in some .sty file and just use it in your preamble.

Q: I want to use

maxima --batch=file

but I need to check automatically that no error occurred. Now unfortunately Maxima does not use the return value (which is always 0, apparently completely un-bothered about the execution of File)—I consider this as a bug?

A: Maxima was not designed for this use, and is hard to use this way because

  1. Maxima often asks questions in the middle of a calculation (e.g. Is x>0?).
  2. Maxima errors are not designed for consumption by other programs.
  3. Some errors are not clean.
  4. Startup time is long.

We’re aware of these issues, and would like to fix them, but they are not fixed yet. You are welcome to work on them with us.

Q: Whenever I finish using Maxima, I find many files cluttering up my $HOME directory, such as maxout.gnuplot and maxout.mgnuplot. How can I configure Maxima/GNUplot to store it’s temp files in ~/.maxima or /tmp?

A: Output files for plotting are supposed to be written to the directory named by maxima_tempdir.

Maxima is supposed to read $HOME/.maxima/maxima-init.mac when it is launched. Try maxima_tempdir:"/tmp" or whatever.

Syntax and Programming Issues

Q: Maxima reports a Lisp error wnen I’m trying to use the symbol called lambda.

A: If you want to see the greek letter lambda, use the symbol %lambda. The symbol lambda has a special meaning to Maxima, related to the lambda calculus.

texput(%lambda, "\\lambda");

tells Maxima to emit \lambda in TeX output.

Q: How to add special characters as valid characters for variable names? For example, I want to make it so that it is valid to begin maxima variable names with “$”.

A: declare("xyz", alphabetic) tells the parser to accept characters x y and z as elements of a symbol. However declare("$", alphabetic) doesn’t seem to have the desired effect, probably because the $ expression terminator is hard wired into the parser.

Perhaps you want to use a backslash.

(%i1) \$foo;
(%o1)                                $foo

Q: I want to define function “gradient”. See space between “grad” and “ient” in output %o3 and the wrong result in %o4.

(%i1) load(vect)$
(%i2) ev(express(grad(x+y+z)),diff);
(%o2)                              [1, 1, 1]
(%i3) gradient(p):=ev(express(grad(p)),diff);
(%o3)              grad ient(p) := ev(express(grad p), diff)
(%i4) gradient(x+y+z);
(%o4)                           grad [1, 1, 1]

A: Maxima’s parser is greedy—foobar is parsed into "foo"(bar) when "foo" is a prefix operator. (Try :lisp $% to see what gradient(p):= ... was parsed into.) Backslash defeats recognition of operators—try \gradient instead. (Likewise a\+b is a symbol, not an expression.)

Q: For some reason, something I did forced me into the Lisp debugger, and now I can’t get out of it: my terminal is full of lines like

Break 5 [27]>

and helpful commands (like :top, :R1) which don’t in fact seem to have any effect at all, except for driving me deeper into the debugger. And even Ctrl-C doesn’t work. How do I get out of this loop and back to the Maxima prompt? I’m using Clisp.

Clisp should have printed out something like:

The following restarts are available:

ABORT          :R1      Abort debug loop
ABORT          :R2      Abort debug loop
ABORT          :R3      Abort debug loop
ABORT          :R4      Abort debug loop
MACSYMA-QUIT   :R5      Maxima top-level

From this, :r5 is the way to get to Maxima top-level, so enter :r5. For me, I get back to Maxima’s repl. The prompt “Break 5” says you’re 5 levels deep in the debugger, so :r1 just brings you one level up. The prompt should change to “Break 4”.

Q:

if foo then do bar

produces an infinite loop.

A: do is the loop operator in Maxima. do bar is equivalent to while true do bar.

You should write simply

if foo then bar

Q: How do convert a list to a sum, like from [a,b,c] to a+b+c?

A:

(%i1) apply("+",[a,b,c]);
(%o1) c+b+a

N.B. To name an operator with syntax, whether infix, prefix, or otherwise, quote it.

Instead of apply, you can use xreduce, rreduce, lreduce, or tree_reduce (see user documentation).

Q: Are threre any way to insert elements in the middle of a list?

insert (e,l,p) := (
  if p=1 or l=[]
    then cons(e,l)
  else
    cons(first(l),insert(e,rest(l),p-1))
)$

N.B. Maxima does not provide functions for in-place modification of lists because lists are considered as values, not as objects.

Q: What is the command to remove only the first element from a list?

A: Maxima has no destructive list operations—it mostly treats lists as immutable—, so you cannot remove an element from a list.

(%i1) a:[1,2,3]$
(%i2) delete(1,a);
(%o2)                               [2, 3]
(%i3) rest(a);
(%o3)                               [2, 3]
(%i4) a;   /* list is unmodified */
(%o4)                              [1, 2, 3]

What you can do is take the rest of a list and assign it.

(%i5) a: rest(a);
(%o5)                               [2, 3]

This does not modify other variables looking at the same list value:

(%i6) a:[1,2,3]$
(%i7) b:a$
(%i8) a: rest(a)$
(%i9) a;
(%o9)                              [2, 3]
(%i10) b;
(%o10)                             [1, 2, 3]

Q: How to add my own Lisp-function to Maxima?

A: You can do something like (defun $my_function (x) ...). Then in Maxima my_function(x) will work.

Or you could also do ?my_function(x). The ? before my_function means to call the Lisp function that is literally named my_function.

For the return value, you need to make your Lisp function return forms in the right way. Perhaps the simplest way is to type in the desired expression to Maxima. Then use :lisp $%onnn to look at how Maxima represents it. So x+1 is ((mplus) $x 1).

Q: Is there something like pattern matching implemented in Maxima?

A: Yes. The stuff for rules implements some pattern matching.

Also op can help you with identifying the type.

(%i1) op(a + b + c);
(%o1)                                 +

See the documentation for op, defrule, tellsimp, and related functions.

Q:

(%i1) p: 2*x*y;
(%o1)                              2 x y
(%i2) length(p);
(%o2)                                3
(%i3) q: -2*x*y;
(%o3)                             - 2 x y
(%i4) length(q);
(%o4)                                1

Is this a bug?

A: Not a bug, intentional behavior that Maxima goes out of its way to provide.

Maxima treats 2*x*y as "*"(2,x,y) and -2*x*y as "-"("*"(2,x,y)).

If you prefer, you can use the internal representation, in which the second one is indeed "*"(-2,x,y):

(%i1) length(-2*x*y),inflag:true;
(%o1)                                 3

Maxima uses the external representation for all functions which examine the form of an expression, including length, args, op, part, map, etc. To use the internal form, bind inflag to true.

To show the structure of expressions, try:

struct(ex):=
  if errcatch(op(ex))=[]
    then ex
  else funmake(concat("<",op(ex),">"),map(struct,args(ex)))$

Normal user-visible structure

map(struct,[1,-1,x,1/x,-1/x,2*x,-2*x,x-1]) =>
  [1, <->(1), x, <//>(1, x), <->(<//>(1, x)), <*>(2, x), <->(<*>(2, x)), <+>(x, <->(1))]

Internal structure

map(struct,[1,-1,x,1/x,-1/x,2*x,-2*x,x-1]),inflag:true =>
  [1, - 1, x, <^>(x, - 1), <*>(- 1, <^>(x, - 1)), <*>(2, x), <*>(- 2, x), <+>(- 1, x)]

Q: Please suggest how I may sort columns of a matrix by some column values.

A: How to sort rows:

p: matrix([1,2,3],[1,3,2],[3,1,2])
sortbyrow2(p) :=
  apply('matrix,sort(args(p), lambda([a,b], a[2]<b[2])))$

sortbyrow2(p) =>
    matrix([3,1,2],[1,2,3],[1,3,2])

How to sort columns:

transpose(sortbyrow2(transpose(p)))

Q: I want to define a function make_fn(n) which should return a function fn(x), where the way in which fn(x) is computed depends on the initially provided n.

I tried it as follows:

make_fn(n) := lambda([x],n+x)$

A: Try using funmake:

(%i1) make_fn(n) := funmake('lambda, [[x], n+x])$
(%i2) make_fn(3)(w);
(%o2) w+3

There is also variants that use subst or buildq.

Though second-order functions are certainly a powerful and appropriate technique for many applications, often in Maxima it is easier to use expressions than functions.

Miscellaneous

Q: Is there any printed textbook about Maxima?

A: There is “The Maxima Book” which is linked from the Maxima documentation web page, but I don’t recommend it. Maybe someone will write another book.

Q: Does anybody know if Maxima can be implemented on Pocket PC devices?

A: Maxima could be installed on many Linux devices. As for Microsoft Windows platforms, I know only about the port by Rainer Keuchel: http://www.rainer-keuchel.de/wince/maxima-ce.html but it’s very old and not maintained.

Q: the manual seems to encourage use of foobar(…) but version x.y.z does not know about this function!

A: Looks like you have to enter load(foobar) first.

The reference manual ought to mention that.

Q: How to build Maxima via a Lisp-only process?

A: There is a file INSTALL.lisp which describes how to build Maxima via a Lisp-only process (based on mk defsystem or whatever it’s called).

INSTALL.lisp should be in the distribution tarball, or, failing that, it’s in the top-level Maxima directory in CVS.

The Lisp-only build doesn’t build XMaxima and I don’t know if it arranges things appropriately to enable plotting via Gnuplot.

Q: Is there a way to use matchdeclare() to require that a pattern parameter is both an integer and greater than some value?

(%i1) matchdeclare (q, integerp)$
(%i2) matchdeclare (q, is (q>4))$

A: Only the last matchdeclare is significant—any previous declarations are thrown away. Also, the predicate has to be the name of a function, a lambda expression, or a function call missing its final argument e.g. freeof(x) which is called as freeof(x, some_expr) when the rule is attempted. matchdeclare doesn’t automatically construct a function from a general expression such as is(m > 4).

(%i1) pred1(x) := is (integerp(x) and x>4)$
(%i2) matchdeclare (q,pred1)$

Q: I am using the vi 7.x as a text editor for my source code. I give all my Maxima scripts a .mac extension. However, I have the impression that the syntax highlighting is not really good or working…

A: There is a maxima.vim syntax colorization file in vim 7.

Most of the scripts in the Maxima distribution are named .mac. A few are named .dem (demonstration). Here is my ~/.vim/filetype.vim which tells vim to colorize those files:

if exists("did_load_filetypes")
  finish
endif
augroup filetypedetect
  au! BufRead,BufNewFile *.mac      setfiletype maxima
  au! BufRead,BufNewFile *.dem      setfiletype maxima
augroup END

Vim doesn’t associate .mac with Maxima by default because .mac was already used by some assembler language (if I recall correctly).

See also Maxima syntax highlighting for Kate and KWrite.