0

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?

KevinKim
  • 1,382
  • 3
  • 18
  • 34
  • I'm not entirely getting your problem. What is the `F`-call and I suppose, your code is `[x,fval] = ...`, which you call in a file named `Test.m`, right? Can you provide the values of `lambda` and `mu` for us? – max Mar 26 '20 at 09:48
  • @max my bad, sorry. lambda = mu =1. Yes, you are right. My code is the `[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 it generates error. I will edit the post. Thanks! – KevinKim Mar 26 '20 at 09:57
  • `u` also misses (sorry, didn't saw this before) – max Mar 26 '20 at 10:27
  • At `t = 10`, the derivative `F'(10)`is almost zero, so `fzero` do not know in which direction it should go to reach the root. With a better starting value (`x0`=0.1), `fzero` output `0.41554` on my session. – obchardon Mar 26 '20 at 13:51
  • @obchardon So, this `fzero` highly depends on the initial starting point? even if theoretically one can show that the root is unique? I have also tried to square my entire function `F` and use `fmincon` on `F^2`. That one seems less dependent on the initial starting point. – KevinKim Mar 26 '20 at 13:54
  • `fzero` is a numerical method, a numerical method as nothing to do with analysis so `fzero` do not know that there is only one root, so yes `fzero` highly depends on the initial starting point. – obchardon Mar 26 '20 at 14:00
  • Noticed that `fmincon` also use a numerical method so it also highly depends on the initial starting point. Finding the roots of a function is not an easy task for a computer. – obchardon Mar 26 '20 at 14:15

1 Answers1

1

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)

max
  • 3,915
  • 2
  • 9
  • 25
  • `normcdf(-arg(t))-u`, this part is wrong I am afraid...the things inside `normcdf` are not different by a sign. You can copy and paste my "brutal force" code, and the curve does cross 0. – KevinKim Mar 26 '20 at 12:55
  • @ftxx would you then please correct your example, as this is coming from your code: `@(t) normcdf(sqrt(lambda/t)*(t/mu-1))+...*normcdf(-sqrt(lambda/t)*(t/mu+1))..` – max Mar 26 '20 at 13:27
  • @ftxx BTW, the answer still holds. The function is not defined for negative values of `t` as the `art()` becomes complex. I will updated the figure as soon as you update your example code – max Mar 26 '20 at 13:32
  • Inside the 2nd `normcdf` it is `-sqrt(lambda./t).*(t/mu+1)` but in your code, it is `-sqrt(lambda./t).*(t./mu-1)`. So `+1` vs. `-1`. Am I missing sth? The 1st part inside the `normcdf`, we hare the same – KevinKim Mar 26 '20 at 13:48