当前位置: 首页 > 程序语言 >求助一个数值积分问题,用matlab的quadgk函数来计算,谢谢!

求助一个数值积分问题,用matlab的quadgk函数来计算,谢谢!

作者 wgdd
来源: 小木虫 1750 35 举报帖子
+关注

现求解一个半无限震荡积分问题:

被积函数为s*besselj(0,s*R)./(s.^2-k0^2), 积分区间为[0,inf], 其中R=6, k0=10, 均为已知常数。我打算用两种方法来计算这个问题。(解析解和数值解)

方法一,解析解
对于这样一个积分,可以查到数学手册中积分结果为pi*i/2*besselh(0,1,k0*R). 经过计算可以得到结果为-0.0744。
计算程序
% analytic method
k0=10;
R=6;
p_inc=real(pi*1i/2*besselh(0,1,k0*R));

方法二,采用matlab中的quadgk函数进行计算,由于k0为奇点,所以将积分区间分为两块。[0,k0]和[k0,inf],计算结果为0.0390
计算程序
% numerical method
k0=10;
R=6;
p_f=@(s)(s).*besselj(0,s*R)./(s.^2-k0^2);

p1=quadgk(p_f,0,k0);
p2=quadgk(p_f,k0,inf);

p_incN=p1+p2;

求各位大牛支招,为啥两种方法差别这么大呀~正确的数值积分方法应该是怎么样的呢?十分感谢! 返回小木虫查看更多

今日热帖
  • 精华评论
  • Mr__Right

    如果你查到的解析正确,结果应该是复数-0.07439126817-0.1436835739*I,
    所说的 数学手册 具体是哪本?

    此外,你是只要实部吗?

    积分如果得到的是复数,……一般的应用情况下没有继续纠结下去的必要了

    这个例子,用matlab里面默认的方法似乎是不行的;
    因为被积函数存在s=10的间断点;

    对Bessel函数进行积分本来就波动不容易,加上存在被积函数为正负无穷大的间断点,又是在包含无穷大的区间上积分,
    这个问题的难度系数:非常高。

  • wgdd

    引用回帖:
    2楼: Originally posted by Mr__Right at 2017-07-30 07:19:33
    如果你查到的解析正确,结果应该是复数-0.07439126817-0.1436835739*I,
    所说的 数学手册 具体是哪本?

    此外,你是只要实部吗?

    积分如果得到的是复数,……一般的应用情况下没有继续纠结下去的必要了

    这 ...

    谢谢您的回复。

    这个解析解的公式我是在一本声学手册的附录中找到的。
    是这本书:https://www.win.tue.nl/~sjoerdr/papers/boek.pdf   附录中公式(D.69)。
    关于这类公式更详尽的证明与推广,也有一篇论文进行了说明:http://aip.scitation.org/doi/full/10.1063/1.3596359

    恩,我只要实部,因为算出来的是声压嘛,取实部才有意义~

    您说的积分得到复数,就没有继续纠结下去的必要了,,这个我不太理解哎?是算出来的是复数就不准吗?

    因为考虑到s=10的奇异性,我把积分区间分成了两段[0,k0]和[k0,inf],我查了下,matlab是推荐采用quadgk时这么处理奇异积分的。

    恩恩,Bessel函数在快速震荡,这个积分确实很难求解。。。不知道您有何高见呢

  • Mr__Right

    引用回帖:
    3楼: Originally posted by wgdd at 2017-07-30 08:58:37
    谢谢您的回复。

    这个解析解的公式我是在一本声学手册的附录中找到的。
    是这本书:https://www.win.tue.nl/~sjoerdr/papers/boek.pdf   附录中公式(D.69)。
    关于这类公式更详尽的证明与推广,也有一篇论文进行 ...

    忙别的事情,没太多时间。

    从文献给出的情况看,如果 k 的虚部为0 (公式中没有给出任何结果)

    很可能这个积分是不存在的 (或者发散的)

  • wgdd

    引用回帖:
    4楼: Originally posted by Mr__Right at 2017-07-30 13:11:22
    忙别的事情,没太多时间。

    从文献给出的情况看,如果 k 的虚部为0 (公式中没有给出任何结果)

    很可能这个积分是不存在的 (或者发散的)...

    仿照我上面贴出来的论文中的证明过程,可以用留数定理证明出来,当k的虚部为0,解析解就是上面我写的那个的吧。
    十分感谢~

  • cooooldog

    引用回帖:
    5楼: Originally posted by wgdd at 2017-07-30 13:35:20
    仿照我上面贴出来的论文中的证明过程,可以用留数定理证明出来,当k的虚部为0,解析解就是上面我写的那个的吧。
    十分感谢~...

    可以做出来。不过,用Matlab计算的话,内置的数值积分算法都不好用。
    所以,很多代码都需要改写,工作量比较大。

    考虑到强行数值计算比较繁琐,化繁为简的话,楼主自己的解析式子的结果没什么问题,拿着用就是了

  • wgdd

    引用回帖:
    6楼: Originally posted by cooooldog at 2017-07-30 16:31:37
    可以做出来。不过,用Matlab计算的话,内置的数值积分算法都不好用。
    所以,很多代码都需要改写,工作量比较大。

    考虑到强行数值计算比较繁琐,化繁为简的话,楼主自己的解析式子的结果没什么问题,拿着用就是 ...

    大牛您好,谢谢回复。我在研究一个问题,上面那个积分是简化了以后的。所以存在解析解,稍微变一下就没有解析解了。我想着把这个积分的程序搞定了(能跟解析解算得差不多准了)后,用数值解来解决一些复杂的问题~

    如果可以的话,能给些建议吗?应该如何改写,哪儿有类似的程序呢?可以学习一下。

  • cooooldog

    引用回帖:
    7楼: Originally posted by wgdd at 2017-07-30 16:35:37
    大牛您好,谢谢回复。我在研究一个问题,上面那个积分是简化了以后的。所以存在解析解,稍微变一下就没有解析解了。我想着把这个积分的程序搞定了(能跟解析解算得差不多准了)后,用数值解来解决一些复杂的问题~
    ...

    通用的解决方案比较难,大部分只能针对特定的问题找特定的解法。

    如果能够通用,mathworks 大概率会集成到某个现成的函数上

猜你喜欢