# Polynomial Power of Octave

This tenth article of the mathematical journey through open source, explains advanced polynomial mathematics in octave.

<< Ninth Article

## 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.

Eleventh Article >>

Send article as PDF

# Get Set with Polynomials in Octave

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 2x3 + 3x2 + 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: x2 + y2 = 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
x2 + y2 = 25 ⇒ 9x2 + 9y2 = 9*25 ⇒ 9x2 + (24 – 4x)2 = 225 ⇒ 25x2 – 192x + 351 = 0
For L2
x2 + y2 = 25 ⇒ x2 + (5√2 – x)2 = 25 ⇒ 2x2 – 10√2x + 25 = 0
For L3
x2 + y2 = 25 ⇒ x2 + (36 – 6x)2 = 25 ⇒ 37x2 – 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 >>

Send article as PDF

# Solve Non-linear Equations using Linear Algebra

This eighth article of the mathematical journey through open source, solves non-linear equations using linear algebra in octave.

<< Seventh Article

Hope you have found out the vegetable prices from the vegetable seller, who had placed various equal priced stacks for sell at ₹30. Recall: One stack had 4 lemons, 7 cucumbers, 9 tomatoes. Another had 2 lemons, 5 cucumbers, 27 tomatoes. And the third had just 9 cucumbers & 15 tomatoes. Prices would be ₹2.00 per lemon, ₹2.50 per cucumber, ₹0.50 per tomato, computed as follows:

\$ octave -qf
octave:1> N = [
> 4 7 9
> 2 5 27
> 0 9 15
> ];
octave:2> inv(N) * [30; 30; 30]
ans =

2.00000
2.50000
0.50000

octave:3>

## Polynomial solving

Note that though in 3 variables, even this was a linear equation. How about solving higher order polynomial equations, meaning of squares, cubes, … of the variables. Say, we want a solution for x in x3 + 3x2 + 3x + 1 = 0. Simple! First define a function for this polynomial. And, then use the function solver fsolve() to solve it, as follows:

\$ octave -qf
octave:1> function y = f(x)
> y = x^3 + 3*x^2 + 3*x + 1;
> endfunction
octave:2> [x, fval, info] = fsolve(@f, 0)
x = -0.99999
fval = 0
info =  1
octave:3>

This indicates the value of x as -0.99999 ≈ -1 as the solution to the function f(x), yielding a function value of 0, with info = 1 indicating that solution is obtained. And you may verify the answer by calling the function f with the variable x as f(x) on the octave prompt. The second parameter in fsolve() is the initial guess of the solution.

## Geometry solving

With the power in hand, why not solve more complex geometric problems? Last time we found the intersection point of two straight lines. How about intersection of a straight line and a circle? Let us have the following straight line and circle, defined in the Cartesian coordinate system, i.e. the x-y system:
4x + 3y = 24
x2 + y2 = 25

To be able to solve it using fsolve(), let’s consider the different variables x & y as fields of a vector X, say x as X(1), y as X(2). Then, the equations can be re-written as follows:
4 * X(1) + 3 * X(2) = 24
X(1)^2 + X(2)^2 = 25
and hence could be solved using fsolve() as follows:

\$ octave -qf
octave:1> function Y = F(X)
> Y(1) = 4 * X(1) + 3 * X(2) - 24;
> Y(2) = X(1)^2 + X(2)^2 - 25;
> endfunction
octave:2> [Y, Fval, info] = fsolve(@F, [0; 0])
warning: matrix singular to machine precision, rcond = 0
warning: attempting to find minimum norm solution
warning: dgelsd: rank deficient 2x2 matrix, rank = 1
Y =

3.0000
4.0000

Fval =

0.0000e+00   2.6691e-07

info =  1
octave:3>

So, (3, 4) is the intersecting point – can be verified by substituting back into the above equations.

## Solve it

Equipped with this knowledge, here’s a couple of teasers for your brain:

1. Find three numbers, product of which is 60; sum of their squares is 50; and their sum is 12.
2. A sage came to a temple with some flowers and dipped all of them into the first magical pond of the temple and got those back, squared. Then, he offered some of those flowers in the temple and dipped the remaining flowers into the second magical pond to get those back, doubled. Then, he again offered the same number of flowers, as offered earlier, and dipped the remaining flowers into the third magical pond to get those back, tripled, which he took back with him as prasadam. Now, the number of flowers he took back with him, is same as in each one of his offerings. Also, what he took back with him is thrice the number of flowers he came with to the temple. How many flowers did he come in with?

If you think, you have got the octave code for solving the above, you may post the solution in the comments below. And as we move on, we would get into specifically playing with polynomials.

Ninth Article >>

Send article as PDF

# Solve Puzzles using Linear Algebra

This seventh article of the mathematical journey through open source, solves puzzles using linear algebra in octave.

<< Sixth Article

Matrix Maths is what is formally called Linear Algebra. We have gone through its basics in the fifth article. Now, we shall apply that to practical usage. What better than solving puzzles using the same.

## Purchase Solving

Shrishti purchased 24 pencils and 12 erasers for ₹96. Divya purchased 20 pencils and 15 erasers for ₹100. What are the prices of the pencil & the eraser?

Assuming that ‘p’ is the price for pencils & ‘e’ is the price for erasers, we have the following two equations:
24 * p + 12 * e = 96
20 * p + 15 * e = 100
Hence, we could get the values of ‘p’ & ‘e’ by solving these equations. Converting them into linear algebra form, they can be re-written using matrix multiplication as:
┏          ┓┏    ┓    ┏        ┓
┃24 12 ┃┃ p ┃    ┃   96 ┃
┃20 15 ┃┃ e ┃ = ┃ 100 ┃
┗          ┛┗    ┛    ┗        ┛
which is Ax = b, ‘x’ being the vector with variables ‘p’ & ‘e’. Hence, we need to solve for x, which is given by: x = A-1b. Using octave:

\$ octave -qf
octave:1> A = [
> 24 12
> 20 15
> ];
octave:2> b = [
> 96
> 100
> ];
octave:3> x = inv(A) * b
x =

2.0000
4.0000

octave:4>

Hence, p = ₹2 and e = ₹4, i.e. each pencil costs ₹2 and an eraser costs ₹4. You may check by putting back these values in our problem statement. Isn’t that cool?

## Geometry Solving

How about finding the intersection point of two straight lines? Let us have the following 2 straight lines, defined in the Cartesian coordinate system, i.e. the x-y system:
4x + 3y = 24
3x + 4y = 25

Similar to the earlier problem, the intersecting point could be obtained as follows:

\$ octave -qf
octave:1> A = [
> 4 3
> 3 4
> ];
octave:2> b = [
> 24
> 25
> ];
octave:3> X = inv(A) * b
X =

3.0000
4.0000

octave:4>

So, (3, 4) is the intersecting point. Want to see it visually. For that, we would just need to rewrite the straight line equations as follows:
y = (24 – 4x) / 3
y = (25 – 3x) / 4
And then here goes the code:

octave:1> x=-10:0.01:10;
octave:2> plot(x, (24 - 4*x)/3, "b.", x, (25 - 3*x)/4, "g.");
octave:3>

Figure 5 shows the plot generated by the above code.

Figure 5: Intersection of straight lines

## Solve it

Equipped with the puzzle solving basics, here’s one for your brain: A vegetable seller has placed various equal priced stacks for sale at ₹30. One stack has 4 lemons, 7 cucumbers, 9 tomatoes. Another has 2 lemons, 5 cucumbers, 27 tomatoes. And the third has just 9 cucumbers & 15 tomatoes. Can you compute the price of each vegetable?

Hint: Assume the price of lemon, cucumber, tomato as ‘l’, ‘c’, ‘t’, and then form the 3 equations in three variables.

If you think, you have got it, you may post the solution in the comments. And as we move on, we will get into some different kind of puzzle solving.

Eighth Article >>

Send article as PDF