Tag Archives: Computer Algebra

Polynomials in Maxima

This sixteenth article of the mathematical journey through open source, demonstrates the polynomial manipulations using Maxima.

<< Fifteenth Article

Polynomials have fascinated mathematicians since ages. Moreover, because of wide variety of their applications, ranging from basic algebra to puzzles to various sciences. Here, we are going to look at some of the polynomial manipulation functions provided by Maxima, and a couple of real world applications using some of them.

Fundamental polynomial operations

Let’s start with a demonstration of the fundamental polynomial operations, like addition, subtraction, multiplication, division. In all these, whenever needed,we’ll use expand() to expand the polynomials, and string() to display the polynomials in a flattened notation.

$ maxima -q
(%i1) p1: x^2 - y^2 + y$
(%i2) p2: -x^2 - y^2 + x$
(%i3) p3: (x + y)^2$
(%i4) string(p1 + p2);
(%o4)                             -2*y^2+y+x
(%i5) string(p1 + p2 + p3);
(%o5)                          (y+x)^2-2*y^2+y+x
(%i6) string(expand(p1 + p2 + p3));
(%o6)                         -y^2+2*x*y+y+x^2+x
(%i7) string(expand(p3 - p1));
(%o7)                            2*y^2+2*x*y-y
(%i8) string(expand(p1 * p2));
(%o8)                   y^4-y^3-x*y^2-x^2*y+x*y-x^4+x^3
(%i9) string(p1 / p2);
(%o9)                     (-y^2+y+x^2)/(-y^2-x^2+x)
(%i10) string(divide(p1, p2));
(%o10)                           [1,y+2*x^2-x]
(%i11) string(divide(x^2 - y^2, x + y));
(%o11)                              [x-y,0]
(%i12) quit();

Note that the / operator just places the polynomials as a fraction, rather then dividing them. And the function divide() actually divides the first polynomial by the second one, yielding a list with the quotient and the remainder of the division. If the division needs to be with respect to a particular variable, that can be passed as the third argument to divide. Check out the variation below, to understand what it means:

$ maxima -q
(%i1) string(divide(x^2 - y^2 + y, x + y, x));
(%o1)                              [x-y,y]
(%i2) string(divide(x^2 - y^2 + y, x + y, y));
(%o2)                            [-y+x+1,-x]
(%i3) quit();

Coefficients of a polynomial

Extracting the coefficients of a polynomial is another basic and common requirement for polynomial manipulations. Maxima provides two slightly different mechanisms. First one, just finds the coefficient of a given variable or its some power, using coeff(). Second one, segregates a polynomial into the coefficient of a given variable or its power, and the remaining terms, using bothcoef().

$ maxima -q
(%i1) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), x^2));
(%o1)                             [2,2*x*y]
(%i2) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), x));
(%o2)                            [2*y,2*x^2]
(%i3) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), y^2));
(%o3)                          [0,2*x*y+2*x^2]
(%i4) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), y));
(%o4)                            [2*x,2*x^2]
(%i5) string(coeff(expand((x + 2 * y)^50), x^20));
(%o5)                   50604606318512743383040*y^30
(%i6) string(coeff(expand((a + b + c + d)^4), a^3));
(%o6)                            4*d+4*c+4*b
(%i7) quit();

Polynomial fractions

Calculating the greatest common divisor (GCD) is one of the very useful operations to simply fractions of polynomials. Other common requirements are extracting the numerator, the denominator, and the highest power. Here goes the function demonstrations for all of these:

$ maxima -q
(%i1) gcd(x^3 + 3*x^2 + 3*x + 1, x^2 + 3*x + 2);
(%o1)                              x + 1
(%i2) string(ezgcd(x^3 + 3*x^2 + 3*x + 1, x^2 + 3*x + 2));
(%o2)                       [x+1,x^2+2*x+1,x+2]
(%i3) string(denom((x + 1)^-3 * (1 - x)^2));
(%o3)                             (x+1)^3
(%i4) string(num((x + 1)^-3 * (1 - x)^2));
(%o4)                             (1-x)^2
(%i5) hipow(expand((x + 1)^3 + (1 - x)^3), x);
(%o5)                                2
(%i6) quit();

Note that the ezgcd() function lists out the remainder polynomials, along with the GCD.

Polynomial fractions could be differentiated using the powerful ratdiff().

$ maxima -q
(%i1) string(ratdiff((x + 1)^-1 * (1 - x)^2, x));
(%o1)                     (x^2+2*x-3)/(x^2+2*x+1)
(%i2) string(ratdiff(1 / (x + 1), x));
(%o2)                         -1/(x^2+2*x+1)
(%i3) string(ratdiff((x^2 - 1) / (x + 1), x));
(%o3)                                1
(%i4) quit();

And, ratsubst() is a powerful expression substitution function, with intelligence. It would dig into the expression to simplify complicated expressions, including trigonometric ones. Check out the %i5, for one of its powerful application. ratsubst(<new>, <old>, <expr>) replaces the <old> expression by the <new> expression in the complete expression <expr>.

$ maxima -q
(%i1) string(ratsubst(u, x^2, x^3 + 3*x^2 + 3*x + 1));
(%o1)                          (u+3)*x+3*u+1
(%i2) string(ratsubst(u, x^2, (x+1)^3));
(%o2)                          (u+3)*x+3*u+1
(%i3) string(ratsubst(u, x^2, (x+1)^4));
(%o3)                       (4*u+4)*x+u^2+6*u+1
(%i4) string(ratsubst(u, x - 1, x^4 - 2*x^2 + 1));
(%o4)                         u^4+4*u^3+4*u^2
(%i5) string(ratsubst(sin(x)^2, 1 - cos(x)^2, cos(x)^4 - 2*cos(x)^2 + 1));
(%o5)                            sin(x)^4
(%i5) quit();

Variable eliminations & equation solving

Very often, we come across N sets of equations in M sets of unknowns, where M >= N. If M = N, then most likely there exists a unique solution. However, if M > N, then there may be many possible solutions, but with some constraints. And in such case, it would be helpful to deduce some simpler set of expressions. This can be achieved using eliminate() of Maxima. Let’s have two polynomial expressions in 3 variables x, y, z to demonstrate the same.

$ maxima -q
(%i1) p1: x^2 + 2*x*y + z^2$
(%i2) p2: x + y + z$
(%i3) string(eliminate([p1, p2], [x]));
(%o3)                            [2*z^2-y^2]
(%i4) string(eliminate([p1, p2], [y]));
(%o4)                         [-z^2+2*x*z+x^2]
(%i5) string(eliminate([p1, p2], [z]));
(%o5)                         [y^2+4*x*y+2*x^2]
(%i6) quit();

Note that in all the above polynomial expressions, the expressions are assumed to be equated to zero. A beautiful application of the above is solving recurrence relations. Let’s say we have a non-linear equation given by Yn+1 = r * Yn * (1 – Yn). And, we want to solve it for a scenario that as n tends to infinity, the value of Y oscillates between two distinct values. This means that Yn+2 = Yn. Hence, we would have two expressions with three unknowns to solve with, namely:

  1. Yn+1 = r * Yn * (1 – Yn)
  2. Yn = r * Yn+1 * (1 – Yn+1)

So, representing Yn+1 with yn1 and Yn by yn, we may use solve() to solve for yn & yn1 in terms of r. But, if rather we want to get a feel of the equation which pops up, in terms of yn, we would have to use eliminate() to eliminate yn1.

$ maxima -q
(%i1) p1: yn1 = r * yn * (1 - yn)$
(%i2) p2: yn = r * yn1 * (1 - yn1)$
(%i3) string(eliminate([p1, p2], [yn1]));
(%o3)          [yn*(r^3*yn^3-2*r^3*yn^2+(r^3+r^2)*yn-r^2+1)]
(%i4) string(solve([p1, p2], [yn, yn1])); /* Just to show the solution */
(%o4) [[yn = (r-1)/r,yn1 = (r-1)/r],
    [yn = -(sqrt(r^2-2*r-3)-r-1)/(2*r),yn1 = (sqrt(r-3)*sqrt(r+1)+r+1)/(2*r)],
    [yn = (sqrt(r^2-2*r-3)+r+1)/(2*r),yn1 = -(sqrt(r-3)*sqrt(r+1)-r-1)/(2*r)],
    [yn = 0,yn1 = 0]]
(%i5) quit()

Seventeenth Article >>

   Send article as PDF   

Expression Simplification using Maxima

This fifteenth article of the mathematical journey through open source, demonstrates the various techniques of expression simplification using Maxima.

<< Fourteenth Article

Expression simplification can be done in a variety of ways. Let’s start with the simple ones, and then move onto more powerful ones.

Real number simplification

A simple simplification is to convert a given number into a rational number using ratsimp(). Below follows the demonstration.

$ maxima -q
(%i1) ratsimp(9);
(%o1)                                  9
(%i2) ratsimp(10.0);
rat: replaced 10.0 by 10/1 = 10.0
(%o2)                                 10
(%i3) string(ratsimp(4.2)); /* string to print it on one line */
rat: replaced 4.2 by 21/5 = 4.2
(%o3)                               21/5
(%i4) string(factor(3.2)); /* string to print it on one line */
rat: replaced 3.2 by 16/5 = 3.2
(%o4)                                2^4/5
(%i5) string(ratsimp(4.3333)); /* string to print it on one line */
rat: replaced 4.3333 by 43333/10000 = 4.3333
(%o5)                            43333/10000
(%i6) string(ratsimp(4.3333333)); /* string to print it on one line */
rat: replaced 4.3333333 by 13/3 = 4.333333333333333
(%o6)                               13/3
(%i7) quit();

Another one is to check whether a number is an integer using askinteger(). And if yes, is it an even or odd number, again using askinteger(). Moreover, asksign() checks for the sign. In case of trying these with unknown variables, Maxima would ask the user the necessary question(s) to deduce the response, and store it for future analysis. For example, askinteger(x) would ask us if x is an integer. And, if we say yes, it can then deduce many more information by itself. Below are some examples:

$ maxima -q
(%i1) askinteger(1);
(%o1)                                 yes
(%i2) askinteger(1.0);
rat: replaced 1.0 by 1/1 = 1.0
(%o2)                                 yes
(%i3) askinteger(1.2);
rat: replaced 1.2 by 6/5 = 1.2
(%o3)                                 no
(%i4) askinteger(-9);
(%o4)                                 yes
(%i5) askinteger(2/3 + 3/4 + 1/6 + 5/12);
(%o5)                                 yes
(%i6) askinteger(-9, even);
(%o6)                                 no
(%i7) askinteger(0, even);
(%o7)                                 yes
(%i8) properties(x);
(%o8)                                []
(%i9) askinteger(x);
Is x an integer?
yes; /* This is our response */
(%o9)                                 yes
(%i10) askinteger(x + 9);
(%o10)                                 yes
(%i11) askinteger(2 * x, even);
(%o11)                                yes
(%i12) askinteger(2 * x + 1, even);
(%o12)                                no
(%i13) askinteger(x, even);
Is x an even number?
n; /* This is our response */
(%o13)                                no
(%i14) askinteger(x, odd);
(%o14)                                yes
(%i15) askinteger(x^3 - (x + 1)^2, even);
(%o15)                                no
(%i16) properties(x);
(%o16)          [database info, kind(x, integer), kind(x, odd)]
(%i17) asksign((-1)^x);
(%o17)                                neg
(%i18) asksign((-1)^(x+1));
(%o18)                                pos
(%i19) asksign((-1)^x+1);
(%o19)                               zero
(%i20) quit();

Maxima simplification uses the concept of properties of symbols. As an example, note the properties of x in the above demonstration at %i8 and %i16.

Complex number simplification

Complex numbers have two common useful forms, namely the exponential and the circular (with sine & cosine) forms. demoivre() converts the exponential to circular form and exponentialize() does the other way round. Using expand() along with them, can simplify complicated looking expressions. Here goes few examples:

$ maxima -q
(%i1) string(demoivre(exp(%i*x^2))); /* %i is sqrt(-1) */
(%o1)                       %i*sin(x^2)+cos(x^2)
(%i2) string(exponentialize(a*(cos(t)))); /* %i is sqrt(-1) */
(%o2)                    a*(%e^(%i*t)+%e^-(%i*t))/2
(%i3) string(expand(exponentialize(a*(cos(t)+%i*sin(t))))); /* %i is sqrt(-1) */
(%o3)                            a*%e^(%i*t)
(%i4) quit();

Expansions and Reductions

As already seen in the previous article, expand() expands an expression completely, by default. However, it can be controlled by specifying the maximum power to which to expand for both the numerator and the denominator, respectively. Using factor() can compact the expanded expressions. Moreover, in many cases, we would like to expand it only with respect to only some variable(s). Say, (x + a + b)^2 should be expanded with respect to x. expandwrt() is meant exactly for that. One example of each of these is shown below.

$ maxima -q
(%i1) string(expand(((x+1)^2-x^2)/(x+1)^2, 2, 0));
(%o1)                      2*x/(x+1)^2+1/(x+1)^2
(%i2) string(factor(((x+1)^2-x^2)/(x+1)^2));
(%o2)                           (2*x+1)/(x+1)^2
(%i3) string(expandwrt((x+a+b)^2,x));
(%o3)                      x^2+2*(b+a)*x+(b+a)^2
(%i4) quit();

Expressions containing logs, exponentials, radicals (powers) can be simplified using radcan(). Rule based simplifications can be achieved using sequential comparative simplification function scsim(). Both of these call for a few examples.

$ maxima -q
(%i1) string(radcan(exp(5 * log(x) + log(3 * exp(log(y) / 4)))));
(%o1)                          3*x^5*y^(1/4)
(%i2) radcan((log(2*x+2*x^2)-log(x))/(log(1+1/x)+log(2*x)));
(%o2)                                1
(%i3) expr: a^2 + b^2 + c^2$
(%i4) eq1: a^2 + 2*a*b + b^2 = 4$
(%i5) eq2: a * b = 6$
(%i6) string(scsimp(expr, eq1, eq2));
(%o6)                              c^2-8
(%i7) quit();

Unlike Octave, Maxima by default doesn’t evaluate its expressions, it only simplifies. What it means is that expressions with integers like cos(1), exp(2), sqrt(3), etc. may remain as is in the most simplified form, instead of evaluating to their respective float numerical values. In such cases, we may force the evaluation by passing the option numer. Similar evaluation can be achieved for predicates, using pred.

$ maxima -q
(%i1) cos(1);
(%o1)                             cos(1)
(%i2) cos(1), numer;
(%o2)                        .5403023058681398
(%i3) sqrt(7);
(%o3)                             sqrt(7)
(%i4) sqrt(7), numer;
(%o4)                        2.645751311064591
(%i4) 1 + 2 > 9;
(%o4)                              3 > 9
(%i5) 1 + 2 > 9, pred;
(%o5)                              false
(%i6) string(%e^%pi < %pi^%e); /* string to print it on one line */
(%o6)                         %e^%pi < %pi^%e
(%i7) %e^%pi < %pi^%e, pred;
(%o7)                              false
(%i8) quit();

Summation Simplifications

Symbolical representation and manipulations of summations can be beautifully achieved using sum() and the option simpsum. Check out the code execution below:

$ maxima -q
(%i1) sum(a * i, i, 1, 5);
(%o1)                                15 a
(%i2) sum(a, i, 1, 5);
(%o2)                                 5 a
(%i3) sum(a * i, i, 1, n);
                                      n
                                     ====
                                     \
(%o3)                              a  >    i
                                     /
                                     ====
                                     i = 1
(%i4) simpsum: true$ /* Enable simplification of summations */
(%i5) string(sum(a * i, i, 1, n)); /* string to print it on one line */
(%o5)                            a*(n^2+n)/2
(%i6) string(sum(i^2 - i, i, 1, n) + sum(j, j, 1, n));
(%o6)                         (2*n^3+3*n^2+n)/6
(%i7) string(factor(sum(i^2 - i, i, 1, n) + sum(j, j, 1, n)));
(%o7)                         n*(n+1)*(2*n+1)/6
(%i8) quit();

Sixteenth Article >>

   Send article as PDF   

Doing Algebra with Maxima

This fourteenth article of the mathematical journey through open source, introduces to doing algebra using Maxima.

<< Thirteenth Article

Maxima is a computer algebra system derived from the earlier Macsyma, which in turn, was written in Lisp.

Algebraic capability of Maxima is the one, which puts it apart from the others. We have worked through basic shell maths, bench calculator bc, and finally the ultra capable octave. But none of them could do algebra with symbols, also referred as symbolic or analytic mathematics. Let’s walk through a few examples to demonstrate what it means.

Getting started

Maxima command line can be started on the shell by typing maxima. -q option to it suppresses the initial welcome message. Typing quit(); on the Maxima command line quits Maxima. Typical Maxima commands shall be terminated either with a semicolon (;) or a dollar ($) – the later suppressing any print from the command. Outputs are numbered and prefaced by %o. Inputs are also numbered, and prefaced by %i. Here are a few factorization examples:

$ maxima -q
(%i1) factor(x^3 - 1);
                                       2
(%o1)                        (x - 1) (x  + x + 1)
(%i2) factor((x + 1)^4);
                                     4
(%o2)                              (x + 1)
(%i3) expand(%o2);
                           4      3      2
(%o3)                     x  + 4 x  + 6 x  + 4 x + 1
(%i4) factor(6!);
                                     4  2
(%o4)                               2  3  5
(%i5) expand(6!);
(%o5)                                 720
(%i6) quit();

factor() and expand() are two complementary functions doing factorization & simplification, respectively. Note the way the powers in the output are printed. For example, 6! is factorized as 24 * 32 * 5. In case, we want the output on a single line, we may use the command string(), which we shall use now onwards. Also, note the usage of output labels as inputs to the commands, namely expand(%o2) in the above example.

Integration and differentiation

Octave has enabled us to do integration and differentiation, but more of definite integrals and differences. Yes, it has indefinite operations for polynomials. But with Maxima, we can just plug in the expressions, as we usually do on our maths notebooks. This is what it means:

$ maxima -q
(%i1) integrate(cos(x),x);
(%o1)                               sin(x)
(%i2) diff(cos(x),x);
(%o2)                              - sin(x)
(%i3) string(integrate(1/(a^2 + x^2)^(3/2),x));
(%o3)                        x/(a^2*sqrt(x^2+a^2))
(%i4) string(diff(log(x), x));
(%o4)                                 1/x
(%i5) y: (x+1)^2$ /* $ suppresses the output */
(%i6) string(diff(y, x));
(%o6)                              2*(x+1)
(%i7) string(integrate(y, x));
(%o7)                            x^3/3+x^2+x
(%i8) quit();

Variable assignments can be done in Maxima using colon (:), as shown in the input %i5, above. There onwards, the variable can be used instead of the original expression, as shown in the inputs %i6 & %i7. Note that the variable assignment assigns the complete symbolic expression, unlike the usual value assignments.

Solving equations

Whether it be linear set of equations, polynomial equations, or for that matter non-linear equations, Maxima can solve most of them with quite ease, and that also in a analytical way. Here are a few common examples: 1) The two solutions of a quadratic equation; 2) Two simultaneous equations in two variables; 3) A non-linear equation with irrational number e. Maxima code demonstration follows:

$ maxima -q
(%i1) string(solve(a*x^2 + b*x + c = 0, x));
(%o1)   [x = -(sqrt(b^2-4*a*c)+b)/(2*a),x = (sqrt(b^2-4*a*c)-b)/(2*a)]
(%i2) string(solve([a*x + b*y = c,d*x + e*y = f],[x, y]));
(%o2)        [[x = -(c*e-b*f)/(b*d-a*e),y = (c*d-a*f)/(b*d-a*e)]]
(%i3) string(solve([exp(x) + 1/exp(x) = a],[x]));
(%o3)      [x = log(a/2-sqrt(a^2-4)/2),x = log(sqrt(a^2-4)/2+a/2)]
(%i8) quit();

Equation solving has been one of the most fascinating as well as challenging tasks for solving mathematical problems, since ages. So, when the computer is there to aid it, we try to go beyond the human limitations. We just saw the well-known formula for the two solutions of a quadratic equation. How about the same for a cubic equation? It looks really complicated and hence we typically don’t use it. But, still you want to see it. Here it goes:

$ maxima -q
(%i1) string(solve(a*x^3 + b*x^2 + c*x + d = 0, x));
(%o1) [x = (-sqrt(3)*%i/2-1/2)*(sqrt(27*a^2*d^2+(4*b^3-18*a*b*c)*d+4*a*c^3-b^2*c^2)/
(2*3^(3/2)*a^2)-(27*a^2*d-9*a*b*c+2*b^3)/(54*a^3))^(1/3)-(sqrt(3)*%i/2-1/2)*(3*a*c-
b^2)/(9*a^2*(sqrt(27*a^2*d^2+(4*b^3-18*a*b*c)*d+4*a*c^3-b^2*c^2)/(2*3^(3/2)*a^2)-
(27*a^2*d-9*a*b*c+2*b^3)/(54*a^3))^(1/3))-b/(3*a),x = (sqrt(3)*%i/2-1/2)*(sqrt(27*
a^2*d^2+(4*b^3-18*a*b*c)*d+4*a*c^3-b^2*c^2)/(2*3^(3/2)*a^2)-(27*a^2*d-9*a*b*c+2*b^3)/
(54*a^3))^(1/3)-(-sqrt(3)*%i/2-1/2)*(3*a*c-b^2)/(9*a^2*(sqrt(27*a^2*d^2+(4*b^3-18*a*
b*c)*d+4*a*c^3-b^2*c^2)/(2*3^(3/2)*a^2)-(27*a^2*d-9*a*b*c+2*b^3)/(54*a^3))^(1/3))-b/
(3*a),x = (sqrt(27*a^2*d^2+(4*b^3-18*a*b*c)*d+4*a*c^3-b^2*c^2)/(2*3^(3/2)*a^2)-(27*
a^2*d-9*a*b*c+2*b^3)/(54*a^3))^(1/3)-(3*a*c-b^2)/(9*a^2*(sqrt(27*a^2*d^2+(4*b^3-18*
a*b*c)*d+4*a*c^3-b^2*c^2)/(2*3^(3/2)*a^2)-(27*a^2*d-9*a*b*c+2*b^3)/(54*a^3))^(1/3))-
b/(3*a)]
(%i2) quit();

And that’s what is the power of Maxima.

Plotting graphs

Along with the algebraic processing, Maxima also supports graph plotting for the algebraic expressions, though we must specify the intervals to plot it. It can do 2-D as well as 3-D plots. The following code shows an example for each of these. And figures 22 & 23 show the respective plots.

$ maxima -q
(%i1) plot2d([sin(x), atan(x), tanh(x)], [x, -5, 5])$
(%i2) plot3d(sin(abs(x) + abs(y))/(abs(x) + abs(y)),[x, -12, 12], [y, -12, 12])$
(%i3) quit();
Figure 22: Output of plot2d()

Figure 22: Output of plot2d()

Figure 23: Output of plot3d()

Figure 23: Output of plot3d()

After this first glimpse of what Maxima can do, we are now ready to move on to detail out the specific topics, like expression simplification, polynomials, etc.

Fifteenth Article >>

   Send article as PDF