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

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 ^{3} + 3x^{2} + 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 ^{2} + y^{2} = 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^{2} + y^{2} = 25 ⇒ 9x^{2} + 9y^{2} = 9*25 ⇒ 9x^{2} + (24 – 4x)^{2} = 225 ⇒ 25x^{2} – 192x + 351 = 0

For L2

x^{2} + y^{2} = 25 ⇒ x^{2} + (5√2 – x)^{2} = 25 ⇒ 2x^{2} – 10√2x + 25 = 0

For L3

x^{2} + y^{2} = 25 ⇒ x^{2} + (36 – 6x)^{2} = 25 ⇒ 37x^{2} – 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.

Pingback: Solve Non-linear Equations using Linear Algebra | Playing with Systems

Pingback: Polynomial Power of Octave | Playing with Systems

Christosc=[1 0 -i]

polyout(c,’x’)

r1=roots(c)

d=[1 0 0 -i]

polyout(d,’x’)

r2=roots(d)

Anil Kumar PugaliaPost authorExcellent Christos. Finally after 6.5 years of the post someone posted a solution. 🙂