# fsolve found wrong solution for easy equation

12 views (last 30 days)
MartinM on 22 Aug 2019
Edited: J Chen on 23 Aug 2019
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

Matt J on 23 Aug 2019
Edited: Matt J on 23 Aug 2019
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.
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

#### 1 Comment

MartinM on 23 Aug 2019
Thanks Matt. Perfect

Matt J on 22 Aug 2019
Edited: Matt J on 22 Aug 2019
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 on 23 Aug 2019
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 on 23 Aug 2019
+1. fsolve is way overkill for a problem like this.