## fsolve found wrong solution for easy equation

### MartinM (view profile)

on 22 Aug 2019 at 15:40
Latest activity Edited by J Chen

### J Chen (view profile)

on 23 Aug 2019 at 15:22
Accepted Answer by Matt J

### Matt J (view profile)

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)
Why?
Thanks

R2013b

Answer by Matt J

### Matt J (view profile)

on 23 Aug 2019 at 2:39
Edited by Matt J

### Matt J (view profile)

on 23 Aug 2019 at 2:39

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;
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

### MartinM (view profile)

on 23 Aug 2019 at 6:05
Thanks Matt. Perfect

Answer by Matt J

### Matt J (view profile)

on 22 Aug 2019 at 18:04
Edited by Matt J

### Matt J (view profile)

on 22 Aug 2019 at 18:06

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

### J Chen (view profile)

on 23 Aug 2019 at 13:53
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

### Catalytic (view profile)

on 23 Aug 2019 at 14:50
+1. fsolve is way overkill for a problem like this.