fsolve features

[Rolf Fabian]

On this list , you'll find frequent questions related to octave's implementation
of  function solver 'fsolve', e.g. last week: 

(1)  >>address@hidden  asked:
>Is there an Octave equivalent the matlab function 'fzero' which gives the
>root of a nonlinear scalar function f(X)?

(2) >>  address@hidden   replied:
 > No, there  isn't, who knows why  nobody has ever implemented  one using the
 > fsolve function provided with octave.  You can use that for your needs.

(3) >> address@hidden   replied
it would be nice to have more information about what 'fsolve' can do.

(4) >>  address@hidden replied :
I'd say is well documented in the octave manual.  What's wrong with it?
                                                 ^^^^^^^^^^^^^^^^^^

(1)   NO there isn't!
(2)   a call to fsolve  like :  [x,info]=fsolve('f',x0)   is adequate for
function   (vector) Y =  f ( vector X  ) 
      >>> where Y and X have same length <<<<
      It's supposed to handle equally sized problems only!

    So  fsolve is NOT ABLE  to solve  for roots of a 
 function   (scalar) y=  f ( vector X  ) 
     e.g. intended to find extrema  of  f(X)   at vanishing  gradient( f(X) )
  
    ( I'm not using matlab at all, but from the discussions, I suspect, that
      this is exactly what 'fzero'  does ... ?? )

    >>  ... who knows why  nobody has ever implemented  one ....



 I know .. Let me tell you why .. !  Besides the problem described above,
 it's because 'fsolve' fails even for the most simple SCALAR problems if
 complex arithmetic is  (maybe implicitely) invoked!

EXAMPLE 1
(A) initialize by real random scalar
[x,RC]= fsolve(f='x^2 +1',randn(1)); #>V2.0.14
    #NOTE: prior to V2.0.14 , you have to  define  f  explicitely:
    #function y=f(x) y=x^2+1; end
                                    #[x,RC]=fsolve('f',randn(1))
x = -0.0011765
RC = 3
This returns  *always*  FAIL return code  RC=3, and nonsense root x, because 
'fsolve'
is not able to recognize that the two possible solutions are purely imaginary, 
i.e.
x1= +i     ,   x2=  -i

(B) but what happens if we initialize by complex random scalar ? i.e.
[x,RC]= fsolve(f='x^2 +1',randn(1)+randn(1)*i );                #>V2.0.14
warning: implicit conversion from complex scalar to real matrix
x = 0.0013174
RC = 3
we get essentially same result as in (A) !


EXAMPLE 2:
roots of a polynomial p=[1,0,0,1] with real coefficients / order 3
( by 'roots(p)'  |-  [-1, 0.50000 + 0.86603i  , 0.50000 - 0.86603i ]
  we already know the answer )
( A )  fsolve real init
[x,RC]= fsolve(f='x^3 +1',randn(1) );                   #>V2.0.14
x   = -1.0000
RC = 1                                          #OK
( B ) let's try  complex init
[x,RC]= fsolve(f='x^3+1',randn(1)+randn(1)*i );                 #>V2.0.14
        warning: implicit conversion from complex scalar to real matrix
x   = -1
RC = 1
( C )  and an imag init
[x,RC]= fsolve(f='x^3+1',randn(1)*i );                          #>V2.0.14
warning: implicit conversion from complex scalar to real matrix
x = -0.0016000
RC = 3
always to fails and not even returns the present real solution!

THE 2 REMAINING COMPLEX SOLUTIONS ARE NEVER RETURNED


EXAMPLE 3:
roots of a 'generalized  polynomial' with real coefficients and 
with NON-INTEGER exponents  close to example (B))

function y=f(x) y=x^pi+1; end
[x,info]=fsolve('f',randn(1))
warning: implicit conversion from complex scalar to real matrix
x = -1.0331
RC = 1

**********************************************************
* Here we even get SUCCESS return code RC=1, but returned 
* x   is far  from being  a ROOT  of  function  f  !!                           
            
**********************************************************

Check:    f(x)     |-      ans = -1.1546e-14 - 4.7669e-01i
which is obviously different from  zero. The only indication, that something
might be wrong comes from the displayed warning .. but such a warning can
also be displayed if  returned root is OK , as in example 2 (B) ... !

**********************************************************
* You'll never get a useful solution from the current implementation 
* of 'fsolve' for this or similar examples
**********************************************************

Because of all these pitfalls of  'fsolve',  I decided to write my my personal 
solver 'cfsolve' ,
returning  correct complex roots for this (or similar)  problem as:
x1 = -0.98999 + 0.14112i
x2 = -0.98999  - 0.14112i
x3 =  0.54030 + 0.84147i
x4 =  0.54030  - 0.84147i

NOTE   from these examples we also learn, that 
        f='x^3  +1'  : has THREE  different  roots (nothing NEW)
                                1 real and 1 complex conjugated pair
                         whereas
 f='x^pi+1'  :  has FOUR   different roots ( 2 complex conjugated pairs)
                 even thought the max exponent ( i.e. pi ) is much CLOSER to 3 
than to 4

It' s obvious, that the single real root of 1st (standard polynomial) problem 
splits
to a conjugated  pair of complex solutions,  when  the max exponent  changes 
from (integer)  3 -> (real) pi.

#------------------------------------------------------------------------------------------------#
Solving for roots of vector/matrix  functions with keeping them as general as 
possible
is even much... much more sophisticated   because the multiplicity of solutions 
may
grow tremendously, even to Inf.

e.g.  let's look at the examples above  now in the matrix sense  f= 'x^3 
+eye(2)'
        'cfsolve' -algorithm  returns as expected
   5.0000e-01 + 8.6603e-01i   1.4990e-18 - 3.3878e-19i
  -2.1984e-19 - 6.1458e-19i   5.0000e-01 + 8.6603e-01i

        but also
                 0.49782 - 0.32433i   0.42819 + 0.68053i
  -0.42858 + 0.67785i   0.50218 + 0.32433i

         and   0.47550 + 0.56074i  -0.65150 - 1.15092i
   0.14439 - 0.29724i   0.52450 - 0.56074i

         and  0.48680 - 0.41490i   0.48158 - 0.17351i
                 0.13985 - 1.35376i  -0.98680 - 0.45113i

        ..... and so on  , where all matrices satisfy f==0.
             i.e. multiplicity of solutions is INFINITY here.

On the other hand  a similar 2x2-matrix function  
 f= 'x^3+1' being interpreted as 'x^3 +ones(2)' by octave
leads to  a finite set of only 3 matrix roots

x1= -0.62996 + 0.00000i  -0.62996 - 0.00000i
       -0.62996 - 0.00000i  -0.62996 + 0.00000i

x2=  0.31498 - 0.54556i  0.31498 - 0.54556i
        0.31498 - 0.54556i  0.31498 - 0.54556i

x3=  0.31498 + 0.54556i  0.31498 + 0.54556i
        0.31498 + 0.54556i  0.31498 + 0.54556i


CONCLUSION
as long as 'fsolve' fails to be able to return correct complex roots,
or at least a RELIABLE error code, it makes really little sense to
write sophisticated vector/matrix - functions solvers relying on
an unreliable main module ....



---------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.  To ensure
that development continues, see www.che.wisc.edu/octave/giftform.html
Instructions for unsubscribing: www.che.wisc.edu/octave/archive.html
---------------------------------------------------------------------