[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
fsolve features
From: |
Rolf Fabian |
Subject: |
fsolve features |
Date: |
Mon, 5 Jul 1999 12:03:18 +-200 |
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
---------------------------------------------------------------------
- fsolve features,
Rolf Fabian <=