MATLAB Answers

fsolve found wrong solution for easy equation

12 views (last 30 days)
MartinM
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)
Solution for z=0 is not found...
Why?
Thanks

  0 Comments

Sign in to comment.

Accepted Answer

Matt J
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.
  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

  1 Comment

Sign in to comment.

More Answers (1)

Matt J
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

  2 Comments

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

Sign in to comment.

Sign in to answer this question.

Products


Release

R2013b

Translated by