Finding roots of a polynomial

Specific domains Optimization
1

Dear Scilabers

I have tried to solve a polynomial for roots and find the positive (real) root.

In Python the following code works successfully:

import numpy as np
from scipy.optimize import fsolve

c1 = 27.9811
c2 = 30.9583

f = lambda x: x**4 - c1*x**3 + c2*x - 1
result = fsolve(f, 1.0)
print(result)          # This returns 1.0557487789221027

I have tried to replicate this in Scilab as follows:

c1 = 27.9811
c2 = 30.9583

p = %s.^4 - c1*%s.^3 + c2*%s - 1

r = roots(p, 'f')
ir = find(r==abs(r))    // find positive real root, index ir
disp(r(ir))             // This returns 27.981085, not the answer I want

r =

0.0211731 + 0.0243569i
0.0211731 - 0.0243569i
-0.0423195 + 0.i
27.981085 + 0.i

If on the other hand I replace r = roots(p, ‘f’) with:

r = roots(p, 'e')       // find eigenvalues

This returns:
r =

27.941504 + 0.i
-1.0484733 + 0.i
1.0557487 + 0.i
0.032332 + 0.i

One of these solutions is the answer I’m looing for, namely the third line (see comment in Python code).

I just wonder if anyone can explain the differences I am observing here.

The python code contains a starting estimate for fsolve (= 1.0) and finds the nearest root. Is there a similar way in Scilab?

With kind regards,
Claus

2

Yes, use Scilab,'s fsolve:

function y=f(x)
  c1 = 27.9811;
  c2 = 30.9583;
  y=x^4-c1*x^3+c2*x-1;
end

--> r=fsolve(1,f)
 r  = 

   1.0557481
3

Hi Stéphane

Thank you. I didn’t even think to look up if Scilab has a function of the same name.

With kind regards,
Claus