Way to handle singularity,Division by zero

Hi All, I am trying to find out the integral by algorithm but it is singular at right end limit.

U see beta is 0.5511 and at right hand limit 1 the term (x/((alpha*x*x)+1-alpha)).^2becomes 1 and further subtraction by 1 make it zero. I am trying (1-sqrt of eps) but I am not sure how much accurate the answer will be. With quadgk answer is 0.5.

but how sould I handle the singularity. I need an accurate answer to find value of alpha which in turn will be used to create some profiles. Thanks in advance. Amit Kadam

     clc;
     clear;

     syms r;

     Alpha =1.436269843278397;

     H = 1;
     R = 1;
     Theta = 0;

     Beta =o.5516 
     a = Beta;
     b = 1;        %(r./R);

     npoints = 63;
     h = (b - a) / (npoints);

     node_r = sym(zeros(1,npoints));
     for i=1:npoints
     % The number of intervals must be even => number of points must be odd.
     x = a - h + (h.*i);
     if x==1
         x =1-(sqrt(eps));
     end
     func_r= ( ( ( ( (x) ./ ( (Alpha.*(x).*(x) ) + 1 - Alpha ) ).^2 ) - 1 ) .^(-     0.5) );
     node_r(i) = func_r;

     end

     simpson = (h / 3) .* (node_r(1) + node_r(end)) + ((4 * h) ./ 3) .* ..

Answer 1

[Doug Hull]

Could you calculate the denominator, and if the absolute value is within a certain tolerance you branch the code to a default value for the expression?

I am surprised if in floating point you ever get zero, you might get really close.

It is not clear which denominator is going to zero.

[Amit Kadam]

U see beta is 0.5511 and at right hand limit 1 the term (x/((alpha*x*x)+1-alpha)).^2becomes 1 and further subtraction by 1 make it zero. I am trying (1-sqrt of eps) but I am not sure how much accurate the answer will be with quadgk answer is 0.5.