fsolve found wrong solution for easy equation

Dear everyone, I have a stupid problem with fsolve:

my equation is :

wehre a and b are constante. I would like to solve it for z=linspace (0,20,100) (for exemple)

if z=0, the solution is x=1.

i write :

close all
clc
clear all
c=3e8;
r=30e-6;
Aeff=pi.*r.*r;
lambda0=1030e-9;
om0=2.*pi.*c./lambda0;
C0=0.;
eta=4E-1
b2=-4e-28;
b2=b2;
b3=1.5e-40;
g=sign(b2);
T0=100./om0;
Tr=10./om0;
n2=1e-23;
gamma=2.*pi.*n2./(lambda0.*Aeff);

a=8.*b3.*Tr./(15.*T0.^4);
b= 2.*eta.*b3./(3.*gamma.*T0.^2);

sol= fsolve(@(x) x*x-1+a/b*log( (b*x*x -a) / (b-a)  ),0)

return 

%%% I also try with solve 
syms u
eqn = u.*u-1+a./b.*log( (b.*u.*u -a) ./ (b-a)  )==0;% + 2b.*z
solu = solve(eqn,u)

Solution for z=0 is not found...

Why?

Thanks

Answer 1

[Matt J]

Choosing an initial guess of x=0 in fsolve is often a bad idea because the gradient of the objective will frequently be zero there (this is the case for your problem) or contain division-by-zero or log(0) operations.

In the implementation below, I've made a number of improvements including,

1. Choosing an initial guess, x=3, away from zero.
2. Supplying an analytical gradient calculation using SpecifyObjectiveGradient=true
3. Using log1p() to compute your objective with higher precision
4. Scaling the objective by 1e6 to achieve more natural units

and get correspondingly better results.

function solveit

a = 0.049032234129438;
b = 6.209226846510589e-07;

    
    opts=optimoptions('fsolve', 'SpecifyObjectiveGradient',true);
                            
    sol=fsolve(@fun,3,opts)
    
    function [f,df]=fun(x)
        
        f=x^2-1+a/b*log1p( b*(x^2 -1) / (b-a)  ); 
           f=f*1e6;
        
        df=2*x*(1+a/(b*x^2-a));  
           df=df*1e6;
        
    end


end

This produces,

sol =

   1.000000000000032

[MartinM]

Thanks Matt. Perfect

Answer 2

[Matt J]

fzero works better

fun=@(x) x^2-1+a/b*log( (b*x^2 -a) / (b-a) );
[sol,fval]= fzero(fun,0)

sol =

  -1.000000128628277


fval =

     4.488965127631997e-12

[J Chen]

The function is symetric to the y-axis. It has roots at both 1 and -1. The following command finds the root at 1

[sol,val] = fzero(@(x) x*x-1+a/b*log( (b*x*x -a) / (b-a)  ),0.5)

sol =
    1.0000
val =
   4.0381e-12

A post discusses the difference between fzero and fsolve can be found here.

[Catalytic]

+1. fsolve is way overkill for a problem like this.