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 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:
- Yn+1 = r * Yn * (1 – Yn)
- 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()
