MATLAB Answers

Fzero returning wrong root

4 views (last 30 days)
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)));
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
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
hold off


Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 19 Nov 2020
Edited: David Goodmanson on 19 Nov 2020
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?

  1 Comment

Emily Davenport
Emily Davenport on 19 Nov 2020
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.

Sign in to comment.

More Answers (0)


Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!