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

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 *Y _{n+1} = r * Y_{n} * (1 – Y_{n})*. 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

*Y*. Hence, we would have two expressions with three unknowns to solve with, namely:

_{n+2}= Y_{n}- Y
_{n+1}= r * Y_{n}* (1 – Y_{n}) - Y
_{n}= r * Y_{n+1}* (1 – Y_{n+1})

So, representing *Y _{n+1}* with

*yn1*and

*Y*by

_{n}*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()
```

Pingback: Expression Simplification using Maxima | Playing with Systems

Pingback: High School Trigo with Maxima | Playing with Systems

kenvery cool Anil. The recurrence relation you mention is known as the logistics function (https://en.wikipedia.org/wiki/Logistic_map). Maxima is very good at plotting the sequence, for example

y:.1; a:[y]; for i:0 thru 126 do (ev(y:r*y*(1-y),[r=3.4]), a:append(a,[y]));

wxplot2d([discrete,a],[style,lines],grid2d);

Anil Kumar PugaliaPost authorBeautiful. Thanks Ken for sharing.