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,
- Choosing an initial guess, x=3, away from zero.
- Supplying an analytical gradient calculation using SpecifyObjectiveGradient=true
- Using log1p() to compute your objective with higher precision
- 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
