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

2014 11 26
1 个回答
回答已采纳
更新时间：2014 12 4
5 次查看（30 天）

0 个投票

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.

请先登录，再回答此问题。

 采纳的回答

Mike Hosea
Mike Hosea 2014-11-26

1 个投票

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 个评论
显示 4更早的评论 隐藏 4更早的评论

Basit
Basit 2014-11-27
编辑：Basit 2014-11-27
Thanks Mike Hosea. I have posted the complete code of my problem, with values resulting in warning. I thought that the problem is due to arguments of hankel function so i posted only that part of code previously. As you suggested to use integral instead of quadl, i think the argument of integral must be explicitly in kz, it means i have to manually put the all the basic equations in terms of kz in the argument of "integral" command.
Mike Hosea
Mike Hosea 2014-11-27
编辑：Mike Hosea 2014-11-27
Thanks for the additional detail. As for using INTEGRAL, I'm not sure what requirement you're referring to. There's no difference with INTEGRAL in how you call it in this case. All I did was change QUADL to INTEGRAL and it worked fine. INTEGRAL does require a function handle, whereas QUADL would accept a string, but that is the only difference in requirement that I can think of. Perhaps you used INT instead?
1. INTEGRAL seems to have no difficulty, so stop using QUADL and start using INTEGRAL. If your version of MATLAB doesn't have INTEGRAL, use QUADGK. I think the problem with the results lies elsewhere, perhaps with the formulation or scaling. I noticed that for n=2, it looked like -2.59 or so might be a zero, but there is no zero crossing, as it appears to be a root of even multiplicity, which makes it difficult. The real and imaginary parts are vastly different in scale, and it does not appear that fsolve does relative error control, so perhaps it would be best to use fsolve on @(x)real(myfunc(x))*scale and @(x)imag(myfunc(x)) separately, where scale is some constant that scales the real part up into the same range as the imaginary part. Then combine the results in some way.
2. You'll need to use a loop over the various values of n to add them up. But of course the first order of business is getting the right results from FSOLVE.
Basit
Basit 2014-11-27
@Mike Hosea, Thank you very much. The equation was solved for single n values by using QUADGK, as you mentioned. I may need help again after applying the loop.
Basit
Basit 2014-11-29
编辑：Basit 2014-11-29
@Mike Hosea, I have placed loop on n, as you suggested.
- For n=0, the solution is found with ease and seems to be correct, with initial guess (e.g. x0= 0.3)
- For n= -N->N, the initial guess(x0) needs to be modified according to largest harmonic, as we can see, kzn= kz + 2*n*pi/p; so if i put initial guess as x0= 0.3+2*N*pi/p; the equation is solved. I subtract 2*n*pi/p from the solution to get kz for n=0 harmonic. The real part appears to be approximately correct however the imaginary part is incorrect !. I couldn't understand why the eq. is not solved (no solution found, or solver stopped prematurely errors occur) with the same initial guess as in case of n=0 ?
Mike Hosea
Mike Hosea 2014-11-29
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 2014-12-4
编辑：Basit 2014-12-4
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)

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by