Caroli 发表于 2008-8-7 10:55

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

大家好，我的被积函数类似脉冲，在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 发表于 2008-8-7 14:17

回复 楼主 Caroli 的帖子

把极值点去掉

alljoyland 发表于 2008-8-7 16:19

原帖由 Caroli 于 2008-8-7 10:55 发表 http://www.chinavib.com/forum/images/common/back.gif
大家好，我的被积函数类似脉冲，在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 发表于 2008-8-7 19:32

回复 2楼 sigma665 的帖子

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

Caroli 发表于 2008-8-7 19:45

回复 3楼 alljoyland 的帖子

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

sigma665 发表于 2008-8-7 22:10

用eps来代替

Caroli 发表于 2008-8-8 10:45

回复 6楼 sigma665 的帖子

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

sigma665 发表于 2008-8-8 14:18

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

Caroli 发表于 2008-8-8 19:23

回复 8楼 sigma665 的帖子

对不起，还是不明白。您的意思是直接以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 发表于 2008-8-9 16:41

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

Caroli 发表于 2008-8-12 09:09

回复 10楼 sigma665 的帖子

太谢谢你了。
可以顺利算了，可是，结果还是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 发表于 2008-8-12 09:28

回复 11楼 Caroli 的帖子

没有mydbeslh这个函数

Caroli 发表于 2008-8-13 19:12

回复 12楼 sigma665 的帖子

补充两个函数：
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 发表于 2008-8-19 22:25

回复 楼主 Caroli 的帖子

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

Caroli 发表于 2008-8-22 09:52

回复 14楼 alljoyland 的帖子

看到了，可是我没有找个7.4版，只有7.1，所以想试一下在低版本上看看能不能解决？
你知道哪儿能下载到7.4吗?:@)
页: [1]
查看完整版本: 各位高手，用quadl积分时的积分敏感区（或有极点初）该如何处理？