Fzero returning wrong root

2 views (last 30 days)
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

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)

Categories

Find more on Optimization in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!