This tenth article of the mathematical journey through open source, explains advanced polynomial mathematics in octave.
Roots of i
Let’s start with the solution to our previous brain teaser of finding the square and cube roots of the imaginary number i, which boils down to finding all the roots of the equations x2 = i and x3 = i, respectively. In other words, roots of the polynomials x2 – i and x3 – i, respectively, as follows:
$ octave -qf
octave:1> roots([1 0 -i])
ans =
0.70711 + 0.70711i
-0.70711 - 0.70711i
octave:2> roots([1 0 0 -i])
ans =
-0.86603 + 0.50000i
-0.00000 - 1.00000i
0.86603 + 0.50000i
octave:3>
Displaying Polynomial
Before we start exploring polynomials further, here’s a couple of basic octave functions to visualize the polynomials. polyout() displays a polynomial in our usual known format. polyreduce() removes the redundant leading zero coefficients. Check out the following:
$ octave -qf
octave:1> p = [ 0 0 5 3 -9 4 6 ];
octave:2> polyout(p, "x");
0*x^6 + 0*x^5 + 5*x^4 + 3*x^3 - 9*x^2 + 4*x^1 + 6
octave:3> polyreduce(p)
ans =
5 3 -9 4 6
octave:4> polyout(polyreduce(p), "x");
5*x^4 + 3*x^3 - 9*x^2 + 4*x^1 + 6
octave:5>
Polynomial Arithmetic
As octave represents polynomial as vectors, their addition and subtraction are vector addition and subtraction, respectively. However, the polynomial vectors may be of different length. Hence, they need to be made of same length before addition or subtraction. So, let’s write functions to do the complete operations:
function [ q1 q2 ] = equalize(p1, p2)
# Assuming p1 & p2 to be row vectors
m = max(length(p1), length(p2));
q1 = [ zeros(1, m - length(p1)) p1 ];
q2 = [ zeros(1, m - length(p2)) p2 ];
endfunction
function p = polyadd(p1, p2)
[ q1 q2 ] = equalize(p1, p2);
p = polyreduce(q1 + q2);
endfunction
function p = polysub(p1, p2)
[ q1 q2 ] = equalize(p1, p2);
p = polyreduce(q1 - q2);
endfunction
Assuming the above code is put in the file polynomials.oct, the same can be used as follows:
$ octave -qf
octave:1> source("polynomials.oct");
octave:2> polyout(p1 = [ 1 2 1 ], "x");
1*x^2 + 2*x^1 + 1
octave:3> polyout(p2 = [ 1 -1 ], "x");
1*x^1 - 1
octave:4> polyout(polyadd(p1, p2), "x");
1*x^2 + 3*x^1 + 0
octave:5> polyout(polysub(p1, p2), "x");
1*x^2 + 1*x^1 + 2
octave:6>
Interestingly, octave already have the functions for multiplication and division of polynomials, namely conv() and deconv(), respectively. Here’s a demonstration:
$ octave -qf
octave:1> polyout(p1 = [ 1 2 1 ], "x");
1*x^2 + 2*x^1 + 1
octave:2> polyout(p2 = [ 1 -1 ], "x");
1*x^1 - 1
octave:3> polyout(conv(p1, p2), "x");
1*x^3 + 1*x^2 - 1*x^1 - 1
octave:4> polyout(deconv(p1, p2), "x");
1*x^1 + 3
octave:5> [ q r ] = deconv(p1, p2)
q =
1 3
r =
0 0 4
octave:6>
Polynomial Differentiation and Integration
Octave also provides functions for differentiation and integration of polynomials, namely polyder() and polyint(). Here goes an example of differentiation and definite integral using the same:
$ octave -qf
octave:1> polyout(p = [ 1 2 1 ], "x");
1*x^2 + 2*x^1 + 1
octave:2> polyout(polyder(p), "x");
2*x^1 + 2
octave:3> polyout(polyint(p), "x");
0.33333*x^3 + 1*x^2 + 1*x^1 + 0
octave:4> q = polyint(p);
octave:5> polyval(q, 3) - polyval(q, 0) # Definite integral of p from 0 to 3
ans = 21
octave:6>
What’s next?
Given a set of data points a common requirement in fields of physics, statistics, and many others is to fit a polynomial to it. Going further, we’ll explore the power of octave for the same.