fzero, because complex function value encountered during search

Hi my friend, I have one question regarding the fzero. My code is as follows:

p=11;
spdf = @(x) 3.*(50.^3)./(x+50).^(3+1)./(1-(50./(50+50)).^3);
fun2 = @(x,z) min(x,max(0,fzero(@(y) (100-p-x+y).^(-0.5)-2.*z(2)-4.*z(1).*y,0))).^2.*spdf(x); 
fun3 = @(z) integral(@(x) fun2(x,z),0,50,'arrayvalued', true);
fun21 = @(x,z) min(x,max(0,fzero(@(y)(100-p-x+y).^(-0.5)-2.*z(2)-4.*z(1).*y,0))).*spdf(x); 
fun31 = @(z) integral(@(x) fun21(x,z),0,50,'arrayvalued', true);  
fsolve(@(z)[fun31(z)-10; fun3(z)-178.57145], [0.002; 0.04])

This code could run and the answer is:

ans =
   0.000243425504922
   0.046499744202347

But I test the answer. It is not correct. The information before showing the answer is

because complex function value encountered during search.

(Function value at -40.96 is -0.053117-0.73555i.)

Check function or try again with a different starting value.

No solution found.

fsolve stopped because the relative size of the current step is less than the

default value of the step size tolerance squared, but the vector of function values

is not near zero as measured by the default value of the function tolerance.

I know that fzero cannot deal with complex function. But if I changed "fzero" to be "fsolve". The code is really really slow and the error is still "no solution solved". Another way I tried is I put "abs" in my code and the new code becomes:

p=11;
spdf = @(x) 3.*(50.^3)./(x+50).^(3+1)./(1-(50./(50+50)).^3);
fun2 = @(x,z) min(x,max(0,fzero(@(y) abs((100-p-x+y)).^(-0.5)-2.*z(2)-4.*z(1).*y,0))).^2.*spdf(x); 
fun3 = @(z) integral(@(x) fun2(x,z),0,50,'arrayvalued', true);
fun21 = @(x,z) min(x,max(0,fzero(@(y) abs((100-p-x+y)).^(-0.5)-2.*z(2)-4.*z(1).*y,0))).*spdf(x); 
fun31 = @(z) integral(@(x) fun21(x,z),0,50,'arrayvalued', true);  
fsolve(@(z)[fun31(z)-10; fun3(z)-178.57145], [0.002; 0.04])

But the new error is

Equation solved, fsolve stalled.

fsolve stopped because the relative size of the current step is less than the

default value of the step size tolerance squared and the vector of function values

is near zero as measured by the default value of the function tolerance.

Although I know if change it to be abs of my function, it will not make my result correct. But I am still very interested in how to fix the second code as well. Any one has some ideas? Thanks in advance!

Answer 1

[Star Strider]

The ‘new error’ does not appear to be an error. The fsolve function encountered what appears to be a global minimum (or at least a minimum where the values of the functions are alll near zero), making the solution a success.

[James_111]

Thanks, Star! I tested the answer into my equation.

p=11;

spdf = @(x) 3.*(50.^3)./(x+50).^(3+1)./(1-(50./(50+50)).^3); %fun2 = @(x) min(x,max(0,fzero(@(y) (100-11-x+y)-1-2.*y,0))).*spdf(x);

fun1 = @(x) min(x,max(0,fsolve(@(y) abs((100-p-x+y)).^(-0.5)-2.*0.053050508469204-4.*0.000065498629374.*y,0))).*spdf(x);

integral(@(x) fun1(x),0,50,'arrayvalued', true)

It gives the value of 9.674228915861301. It is far away from the target value of 10.

I checked another equation as follows:

p=11;

spdf = @(x) 3.*(50.^3)./(x+50).^(3+1)./(1-(50./(50+50)).^3); %fun2 = @(x) min(x,max(0,fzero(@(y) (100-11-x+y)-1-2.*y,0))).*spdf(x);

fun1 = @(x) min(x,max(0,fsolve(@(y) abs((100-p-x+y)).^(-0.5)-2.*0.053050508469204-4.*0.000065498629374.*y,0))).^2.*spdf(x);

integral(@(x) fun1(x),0,50,'arrayvalued', true)

It gives the value of 173.121815316693. It is also far away from the target value of 178.5715.

So I cannot accept the answer. I understand your point. But I still need to figure it out if there is any way to fix this issue.

[Star Strider]

The fsolve (and other optimization functions) operate within tolerances. Being ‘close enough’ works for them in order that the algorithms stop searching at some point, rather than consuming system resources and entering into a region where floating-point approximation error keeps them from ever converging absolutely on a zero result. See the discussion on Tolerances and Stopping Criteria (link) and Optimization Options Reference (link) to understand them, and change the stopping creiteria you want.

[James_111]

Thanks, Star! I clicked the stopping criteria details. it shows:

fsolve stopped because the relative norm of the current step, 2.070035e-13, is less than

max(options.StepTolerance^2,eps) = 1.000000e-12. The sum of squared function values,

r = 1.183146e-23, is less than sqrt(options.FunctionTolerance) = 1.000000e-03.

Optimization Metric Options

relative norm(step) = 2.07e-13 max(StepTolerance^2,eps) = 1e-12 (default)

r = 1.18e-23 sqrt(FunctionTolerance) = 1.0e-03 (default)

So I thought the default StepTolerance and FunctionTolerance is 1.000000e-6.

I use two forms of codes to change the StepTolerance after checking your website as follows:

options.StepTolerance = 1e-10;

Or

options = optimoptions(options,'StepTolerance',1e-10);

But none of them works. It means I still get the same answer and the same stopiing criteria details. Do you know why?

[Star Strider]

My pleasure.

It appers that fsolve has reached the limit of floating-point precision. It is unlikely to be able to improve on that.

[James_111]

Thanks, Star! I have an idea to solve my issue suddently as my function is monotone. Can you provide me some idea how to realize the following function?

I want to construct a function of x, where z=[z(1) z(2)] is given

For any x,

If (100-10-x+0).^(-0.5)-2.*z(2)-4.*z(1).*0<=0, then fun1(x,z)=0;

elseif (100-10-x+50).^(-0.5)-2.*z(2)-4.*z(1).*50>=0, then fun1(x,z)=50;

elseif fun1(x,z)=fzero(@(y) (100-p-x+y).^(-0.5)-2.*z(2)-4.*z(1).*y,[0,50])

This approach can may make complex function gone. Thank you so much!

[Star Strider]

As always, my pleasure!

I have a problem following the function you described, since I do not understand the relation between ‘x’ and ‘z’.

[James_111]

Thanks. I created a new problem in the following link

https://www.mathworks.com/matlabcentral/answers/443395-if-function-index-exceeds-matrix-dimensions

[Star Strider]

Noted.

Answer 2

[Matt J]

Since you only have two unknowns, you could plot the function (as a surface). Then you can see graphically if a solution even exists, and if so, approximately where your initial guess should be (i.e., where the surface reaches zero).

In any case, I would use fminsearch to solve this problem. fsolve is not a good choice here, because it is not even clear if your function is smooth.

[James_111]

Thanks, Matt! I will look for the difference between fsolve and fminsearch.