Fzero returning wrong root

Hello,

I've tried everything I can think of to get the intersection point of these two curve, however, no matter which method I've used (interp1, fzero, min(abs(P-y)), ...) I get the wrong answer. I can tell because when I plot the curves with the found intersection point it is not the intersection. Can anyone tell what's going wrong for me? Much appreciated.

syms theta
% integrate Fr (no longer a function of r) w/r.t. theta; 
% move constants out of integrand
g(theta) = (cos(theta))^(1/2)*(cos(theta/2))^9;
F_theta = vpaintegral(g,[-pi() pi()],'reltol', 1e-32, 'AbsTol', 0)   % no closed form solution, numerically approximate

% plot P vs Ki

P = @(Ki) (real(1- exp(-Ki*F_theta)));
figure(1)
Ki_interval = [0 10];
fplot(P, Ki_interval, 'b'); title('Probability of Failure vs. Normalized SIF'); 
xlabel('Normalized Ki'); ylabel('Propability of Failure, P');

% Plot Zoomed in P vs. Ki to get an iniital guess for Ki such that P=0.95
figure(2)
fplot(P, Ki_interval, 'b'); 
xlim([1 3]); ylim([0.75,1]);
hold on
fplot(0.95, Ki_interval, 'r')        % graph P=0.95
title('Zoomed in Probability of Failure vs. Normalized SIF');
xlabel('Normalized Ki'); ylabel('Propability of Failure, P');
hold on

Ki_guess = 2;
%y = @(Ki) (0.95);
fun = @(Ki) ((real(1- exp(-Ki*F_theta)))-0.95);
% fplot(fun, Ki_interval,'k')
% hold on
% fplot(0, Ki_interval, 'g') 
intersection_point = fzero(fun, Ki_guess) 

% graph intersection point to check
scatter(intersection_point,0.95,'or','filled')

hold off

Answer 1

[David Goodmanson]

Hello Emily

If you insert this

F_theta = double(F_theta);

just after the calculation of F_theta, you will get the correct result. The natural domain of fzero is doubles, and mixing in vpa caused problems.

I would be remiss if I did not mention that the entire calculation does not have the right feel. You are creating a cdf, which is a real function. Perhaps the calculation is correct and the imaginary part was intentional. But the vast majority of the time doing a calculation that unexpectedly gives a complex result, then just taking the real part and moving on, is incorrect. Is it possible that the integral for F_theta should be from -pi/2 to pi/2?

[Emily Davenport]

Thank you! I did not know I needed the double type for F_theta. Also, yes you are indeed correct I had the wrong bounds! I appreciate your help.