Roots | Playing with Systems

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 x² = i and x³ = i, respectively. In other words, roots of the polynomials x² – i and x³ – 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>

This ninth article of the mathematical journey through open source, deals with polynomial mathematics in octave.

<< Eighth Article

Let’s first solve the earlier puzzles. And then we shall discuss the polynomial power of octave.

Number Puzzle

Find three numbers, product of which is 60; sum of their squares is 50; and their sum is 12. Let the X vector elements X(1), X(2), X(3) be the three numbers. Then, here goes the solution:

$ octave -qf
octave:1> function Y = F(X)
> Y(1) = X(1) * X(2) * X(3) - 60; 
> Y(2) = X(1)^2 + X(2)^2 + X(3)^2 - 50; 
> Y(3) = X(1) + X(2) + X(3) - 12; 
> endfunction
octave:2> [Y, Fval, info] = fsolve(@F, [3; 3; 3]) 
warning: matrix singular to machine precision, rcond = 4.32582e-35
warning: attempting to find minimum norm solution
warning: dgelsd: rank deficient 3x3 matrix, rank = 1 
Y =

   5.0000
   3.0000
   4.0000

Fval =

  -3.2345e-07   1.0351e-07   0.0000e+00

info = 1 
octave:3>

So, the 3 numbers are 5, 3, 4.

Flower Puzzle

A sage came to a temple with some flowers and dipped them into the first pond of the temple to get them squared. Then, he offered some flowers in the temple and dipped the remaining flowers into the second pond to get them doubled. Then, he again offered same number of flowers, as earlier, and dipped the remaining flowers into the third pond to get them tripled and take back with him as prasadam, which was the same number as in each one of his offerings. Now, if he took back thrice the number of flowers he brought. How many did he bring in with him?
Let the x vector elements x(1) and x(2) be respectively, the number of flowers the sage came with and the number of flowers the sage offered each time. So, here goes the solution:

octave:1> function y = f(x)
> y(1) = ((x(1) * x(1) - x(2)) * 2 - x(2)) * 3 - x(2);
> y(2) = x(2) - 3 * x(1);
> endfunction 
octave:2> [x fval info] = fsolve(@f, [10; 10])
x =

    5.0000
   15.0000

fval =

  -2.8791e-06  -1.7764e-15

info =  1
octave:3>

So, the number of flowers the sage came with is 5 and his each offering is of 15 flowers.

Note that in all these solutions the trick is to choose the initial solution close to the original solution, through some approximation work. At times that might be tricky. So, in case we just have polynomial equations and that also in one variable, it can be solved in an easier way, using the polynomial features of octave. In contrast to the earlier method, here we also get all of the multiple solutions for the polynomial.

Playing with Polynomials

Let’s consider the polynomial equation 2x³ + 3x² + 2x + 1 = 0. Then its octave representation and computation of its solutions aka roots would be as follows:

octave:1> P = [2; 3; 2; 1];
octave:2> roots(P)
ans =

  -1.00000 + 0.00000i
  -0.25000 + 0.66144i
  -0.25000 - 0.66144i

octave:3>

So, it being a cubic equation, it has three roots as expected. First one is the real number -1, and the other two are complex conjugates (-1 + sqrt(-7))/4 & (-1 – sqrt(-7))/4. And you may verify the solutions using the function polyval() as follows:

octave:1> P = [2; 3; 2; 1];
octave:2> sols = [-1; (-1 + sqrt(-7)) / 4; (-1 + sqrt(-7)) / 4]
sols =

  -1.00000 + 0.00000i
  -0.25000 + 0.66144i
  -0.25000 + 0.66144i

octave:3> polyval(P, sols)
ans =

   0
   0
   0

octave:4>

This shows that the value of the polynomial P evaluated at each of the 3 solutions is 0. Hence, confirming that they indeed are the solutions.

All set with polynomial basics in octave, let’s solve some puzzles.

Geometry Solving

Last time we found an intersection point of a straight line and a circle. Yes, we just calculated one point – though typically there would be two. It would be one only in case of the straight line being tangent or just touching the circle. And yes it would be zero, if the straight line is not even intersecting it. So now, let’s try these different cases, with the one variable polynomial power.

Let us have the following circle C with radius 5 and centered at origin (0, 0), defined in the Cartesian coordinate system, i.e. the x-y system: x² + y² = 25

And, let us consider the following 3 lines for intersection with the above circle, one by one:

L1: 4x + 3y = 24
L2: x + y = 5√2
L3: 6x + y = 36

To be able to solve for the intersection points of each of these 3 lines with the circle C using roots, the first step is to get polynomials in one variable. For that, we can substitute the value of y in the equation of the circle, in terms of x from each of the line equations, as follows:
For L1
x² + y² = 25 ⇒ 9x² + 9y² = 9*25 ⇒ 9x² + (24 – 4x)² = 225 ⇒ 25x² – 192x + 351 = 0
For L2
x² + y² = 25 ⇒ x² + (5√2 – x)² = 25 ⇒ 2x² – 10√2x + 25 = 0
For L3
x² + y² = 25 ⇒ x² + (36 – 6x)² = 25 ⇒ 37x² – 432x + 1271 = 0

Now, we get the roots of each to get the x co-ordinate of the intersection point.

octave:1> C1 = [25; -192; 351];
octave:2> C2 = [2; -10*sqrt(2); 25];
octave:3> C3 = [37; -432; 1271];
octave:4> roots(C1)
ans =

   4.6800
   3.0000

octave:5> roots(C2)
ans =

   3.5355
   3.5355

octave:6> roots(C3)
ans =

   5.8378 + 0.5206i
   5.8378 - 0.5206i

octave:7>

And the corresponding y co-ordinate could be obtained by substituting the value of x into the corresponding line equations.

For L1, there are 2 different roots 4.68 and 3, implying two intersecting points (4.68, 1.76) and (3, 4).
For L2, there are 2 identical roots of 3.5355 i.e 5/√2, implying just one intersecting point (5/√2, 5/√2).
For L3, the roots are complex, implying that there is no intersecting point in the real world.

Solve it

And finally, here’s one for your brain. Find out the two square roots and the three cube roots of the imaginary number i.

If you think, you have got the octave code for solving the above, post your solution in the comments below. And as we move on, we would have more fun with the polynomials.

Tenth Article >>

Playing with Systems

SysPlay's Blogs

Tag Archives: Roots

Polynomial Power of Octave