How to remove singularity of quadl in which integrand is hankel function

5 views (last 30 days)
Hi, i am using quadl to integrate hankel function in an m file which is further to be passed to fsolve command.
function F=myfunc(kz)
%------Constants-------
L= 4.2;
W=0.396;
b= 5.592;
a=12;
k0= 0.6283;
B1= -2.0726 - 0.0004i;
C1= 0.8285;
E1= 0 + 0.3152i;
G1=-0.2944;
G2= 1.7056;
p= 0.99;
n=2;m=1;
%--- Basic equations---------
kp= sqrt(k0^2-(kz+2*n*pi/p)^2);
kymn1= sqrt(k0^2-(m*pi/a)^2-(kz+2*n*pi/p)^2);
cotterm=(1/kymn1/b).*cot(kymn1.*b);
Sinc_F= (sinc((kz+2*n*pi/p)*W/2))^2;
%----Integral
hank= @(r) besselh(0,2,(kp)*r*L/pi);
int1= @(r) G1*pi*cos(r)-G1.*r.*cos(r)+G2*sin(r);
intanswer= quadl(@(r)(int1(r)).*hank(r),0,pi);
%-----Final Equation ------
F= Sinc_F*(B1*C1*cotterm+E1*intanswer);
In command window i call
>> fsolve(@myfunc,1) % 1 is the initial guess for kz
My problem:
1- for n=0, the equation is solved. For values of n other than 0, e.g. n=2, matlab gives this message. Warning: Maximum function count exceeded; singularity likely. > In quadl at 106
The equation is solved, but the result is not correct (specially the imaginary part of solution which comes to be zero and should not be zero) which i think is due to the above warning.
2- For simplicity i used only single value of n. But actually it should be let's say from -5 -> 5. I have to add up the effect of basic equations and integrals for all values of n, and then put this in final equation to solve. Any suggestions how can i put summation in above equations.
Please help how to resolve this problem. Thanks in anticipation.

Accepted Answer

Mike Hosea
Mike Hosea on 26 Nov 2014
QUADL is obsolete. You haven't supplied any values for k0 and L. Please provide just one value for each of the following inputs so that this does NOT work. Then we can see what the problem is.
k0 = 1;
p = 1;
n = 1;
L = 1;
f = @(kz)integral(@(r)besselh(0,2,(sqrt(k0^2-(kz+2*n*pi/p)^2))*r*L/pi),0,pi);
fsolve(f,1)
  6 Comments
Mike Hosea
Mike Hosea on 29 Nov 2014
I'll have to take a closer look when I get the chance, but zeros can be hard to approximate numerically sometimes. It is a matter of conditioning, and the case I looked at had both an even root and vastly different scaling between the real and imaginary parts. This different scaling is accompanied by lower relative precision in the value of the other part, whether it be real or imaginary. Plus the solver would sometimes quit early (compared to what I wanted) because of the way it was applying the tolerance.
Basit
Basit on 4 Dec 2014
Edited: Basit on 4 Dec 2014
Dear Mike Hosea. Thanks for your comments. The above code has some errors in the constant values. I have rather a general question. Consider n=0 harmonic, the equation of kpn = sqrt(k0^2-kz^2). It is assumed that real(kz)<k0. I gave initial value of kz which is real and less than k0. But when fsolve uses iterations to solve my complete problem (kpn is included in the problem) it puts several complex values of kz for different iterations. Can you suggest some way to put a restriction, so that fsolve tries for only those values of kz for which real(kpn)>0 and also imag(kpn)>0. ? (The resulting solution kz of the problem needs to be a complex number)

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!