CODEWITHSUNDEEP
matlabnumerical-methods
KevinKim
KevinKim

Reputation: 1434

`fzero` function in matlab fails to find the root and keeps on generating error

I am using MATLAB and I want to find the root of an equation F(x)-u=0. Here u=0.2861 and 

F=normcdf(sqrt(lambda/t)*(t/mu-1))+exp(2*lambda/mu)*normcdf(-sqrt(lambda/t)*(t/mu+1)).

The value of lambda and mu are both 1.

I typed the following code

[x,fval] = fzero(@(t) normcdf(sqrt(lambda/t)*(t/mu-1))+exp(2*lambda/mu)*normcdf(-sqrt(lambda/t)*(t/mu+1))-u, 10);

and hope this can help me find the root. I can show mathematically that this equation has unique root. However, I keep on getting the following error

Error using erfc Input must be real and full.

Error in normcdf>localnormcdf (line 128) p(todo) = 0.5 * erfc(-z ./ sqrt(2));

Error in normcdf (line 50) [varargout{1:max(1,nargout)}] = localnormcdf(uflag,x,varargin{:});

Error in Test>@(t)normcdf(sqrt(lambda/t)*(t/mu-1))+exp(2*lambda/mu)*normcdf(-sqrt(lambda/t)*(t/mu+1))-u

Error in fzero (line 363) a = x - dx; fa = FunFcn(a,varargin{:});

Then I did a "brutal force" method.

t = [0:0.001:20];
F = normcdf(sqrt(lambda./t).*(t/mu-1))+exp(2*lambda/mu).*normcdf(-sqrt(lambda./t).*(t/mu+1))-u;
plot(t,F)

I can clearly eyeball that F(t)-u is increasing in t and the root is around 0.4. My question is why fzero does not work in this case and is there a way to make fzero work?

Upvotes: 0

Views: 918

Answers (1)

Max
Max

Reputation: 4045

The problem is that the function does not change sign, which is required as the docs say:

x = fzero(fun,x0) tries to find a point x where fun(x) = 0. This solution is where fun(x) changes sign — fzero cannot find a root of a function such as x^2.

I broke up your code to make it a bit clearer (at least for me). 

lambda = 1;
mu = 1;
u = 1;

% break up function code
arg1 = @(t) +sqrt(lambda./t).*(t./mu-1);
arg2 = @(t) -sqrt(lambda./t).*(t./mu+1);
fnc = @(t) normcdf(arg1(t))+exp(2*lambda/mu).*normcdf(arg2(t))-u;
% call fzero to find the root
% [x,fval] = fzero(fnc, 10);

% plot
x = 0:0.01:10;
plot(x,fnc(x))

The function is not defined for any input t < 0 due to the sqrt in my function handle arg. So if you plot it for values t > 0, you see that it never passes zero.

EDITED: sign mix-up in the arguments. Thx flxx for pointing this out. Plot & code updated. The argument still holds.

fnc(t)

Upvotes: 1

Related Questions