各位高手，用quadl积分时的积分敏感区（或有极点初）该如何处理？

Caroli

大家好，我的被积函数类似脉冲，在0附近很大，然后在0.1以后就会迅速衰减为零。这样用不同积分算下的结果就差异很大：
当pp=quadl('try0806',-0.01,-10)时，结果为：pp =4.963518457903922e+148 -1.527613904435643e+149i，
当pp=quadl('try0806',-0.01,-10)时，结果为：pp =1.701536867603408e+002 -5.236792005319768e+002i
但是，实际计算的物理量不可能应该不是很大，在e2数量级以内吧。
附：程序：
％％％％
function y=try0806(sgm)

kz=9.42477796076938;rr =0.23994497937820;
r0 =0.10000000000000;fy0=0;
r =0.30000000000000;fy=0.78539816339745;
R=0.50000000000000;m=2;
rou=1.22500000000000;c=340;
e0=exp(i*kz*abz);
          jr0=besselj(m,i*kz*r0./tanh(sgm));

          B=-mydbeslh(m,i*kz*R./tanh(sgm)).*(jr0./mydbeslj(m,i*kz*R./tanh(sgm)));

          jr=besselj(m,i*kz*r./tanh(sgm));   

          y=-rou*c/(8*pi)*kz*kz*e0.*B.*jr.*cos(m.*(fy-fy0))./(tanh(sgm).*sinh(sgm))


％％％％％
main0807
pp=quadl('try0806',-0.01,-10)


请指点，到底该如何确定积分区间？？？？谢谢

sigma665

把极值点去掉

alljoyland

原帖由 Caroli 于 2008-8-7 10:55 发表

                               
登录/注册后可看大图

大家好，我的被积函数类似脉冲，在0附近很大，然后在0.1以后就会迅速衰减为零。这样用不同积分算下的结果就差异很大：
当pp=quadl('try0806',-0.01,-10)时，结果为：pp =4.963518457903922e+148 -1.527613904435643 ...

楼主函数中的abz未定义

同时 对于有 奇异的积分 仅限于 端点处,
楼主可以使用 quadgk 函数, 7.4的版本里面有,
低版本不知道是否有,
其次楼主可以参看 奇异积分 数值积分的内容

后来发现其实还有些函数未定义的,
楼主为什么不给至少可以运行的程序呢,

[ 本帖最后由 alljoyland 于 2008-8-7 16:24 编辑 ]

Caroli

嗯
“0”就是极点了，但是极点怎么去啊？就是该怎么去呢？因为去掉（-0.01，0.01）或者（-0.1，0.1）结果就差很多了。请指教！

Caroli

呵呵，不好意思，粘贴的时候把“abz=1.2;”丢掉了。现在应该好了。

sigma665

用eps来代替

Caroli

呵呵，我比较笨，看了半天help也不知道这个eps怎么用到我的程序里。请前辈说地稍微详细一点。万分感谢！

sigma665

你说极值点是0点
把0点扣除，用一极小值eps代替

Caroli

对不起，还是不明白。您的意思是直接以eps（0）代替0吗？因为以前没接触过类似的问题。
本来是pp=quadl('try0806',0,-10)，但是因为“0”是极点，要将“0”抠掉。故要找到一个合适的接近“0”的数，您的意思是接近“0”的数就是“eps”或者“eps（0）”吗？即：
pp=quadl('try0806',eps,-10)或者pp=quadl('try0806',eps（0）,-10)
但是这样老出错。

且：>> eps(0)

ans =

  4.9407e-324
那我直接用‘4.9407e-324’代替”0“，同样也是不行啊
呵呵
麻烦了

sigma665

pp=quadl('try0806',-10,-eps)
程序不完整，没办法调试
而且，你始终没有贴出错误提示
不过，估计应该可以

Caroli

太谢谢你了。
可以顺利算了，可是，结果还是nan。能不能再帮我看看。附：子程序，主程序及运行结果。
子程序：
function y=try0806(sgm)

kz=9.42477796076938;rr =0.23994497937820;
r0 =0.10000000000000;fy0=0;
r =0.30000000000000;fy=0.78539816339745;
R=0.50000000000000;m=2;
rou=1.22500000000000;c=340; abz=1.2;
e0=exp(i*kz*abz);
          jr0=besselj(m,i*kz*r0./tanh(sgm));

          B=-mydbeslh(m,i*kz*R./tanh(sgm)).*(jr0./mydbeslj(m,i*kz*R./tanh(sgm)));

          jr=besselj(m,i*kz*r./tanh(sgm));   

          y=-rou*c/(8*pi)*kz*kz*e0.*B.*jr.*cos(m.*(fy-fy0))./(tanh(sgm).*sinh(sgm))
主程序：
clear all
clc
pp=quadl('try0806',-eps,-10)
运行结果：
Warning: Infinite or Not-a-Number function value encountered.
> In quadl at 101
  In main0807 at 3

pp =

      NaN +    NaNi

sigma665

没有mydbeslh这个函数

Caroli

补充两个函数：
function yy=mydbeslh(m,z)
yy=(besselh(m-1,z)-besselh(m+1,z))/2;


function yy=mydbeslj(m,z)
yy=(besselj(m-1,z)-besselj(m+1,z))/2;
呵呵，谢谢你对我这个小问题的关注。真的，很谢谢你！:handshake

alljoyland

quadgk
我不是提到了这个函数么,似乎用这个就不用管极点了,怎么没看到呢

Caroli

看到了，可是我没有找个7.4版，只有7.1，所以想试一下在低版本上看看能不能解决？
你知道哪儿能下载到7.4吗?:@)