应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » MATLAB/Mathematica » 求助一个数值积分问题，用matlab的quadgk函数来计算，谢谢！
361/4返回列表1234下一页
查看: 3068  |  回复: 35
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

wgdd

木虫 (正式写手)

[求助] 求助一个数值积分问题，用matlab的quadgk函数来计算，谢谢！ 已有1人参与

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

被积函数为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;

求各位大牛支招，为啥两种方法差别这么大呀~正确的数值积分方法应该是怎么样的呢？十分感谢！
回复此楼

» 猜你喜欢
1楼 2017-07-29 19:37:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

Mr__Right

专家顾问 (著名写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
jjdg: 金币+1, 感谢参与 2017-07-30 09:27:48
wgdd(月只蓝代发): 金币+20, 感谢帮助！ 2017-08-03 12:14:30
如果你查到的解析正确，结果应该是复数-0.07439126817-0.1436835739*I,
所说的 数学手册 具体是哪本？

此外，你是只要实部吗？

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

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

对Bessel函数进行积分本来就波动不容易，加上存在被积函数为正负无穷大的间断点，又是在包含无穷大的区间上积分，
这个问题的难度系数：非常高。
赞一下(1人) 回复此楼

» 本帖已获得的红花（最新10朵）

wgdd
文章乃身外之物,要多考虑编辑、审稿人和读者的感受。
2楼2017-07-30 07:19:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

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函数在快速震荡，这个积分确实很难求解。。。不知道您有何高见呢？
赞一下 回复此楼
3楼2017-07-30 08:58:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

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 （公式中没有给出任何结果）

很可能这个积分是不存在的 （或者发散的）
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

wgdd
文章乃身外之物,要多考虑编辑、审稿人和读者的感受。
4楼2017-07-30 13:11:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wgdd

木虫 (正式写手)
送红花一朵
引用回帖:
4楼: Originally posted by Mr__Right at 2017-07-30 13:11:22
忙别的事情，没太多时间。

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

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

仿照我上面贴出来的论文中的证明过程，可以用留数定理证明出来，当k的虚部为0，解析解就是上面我写的那个的吧。
十分感谢~
赞一下 回复此楼
5楼2017-07-30 13:35:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cooooldog

铁杆木虫 (著名写手)

ส็็็

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
jjdg: 金币+1, 感谢参与 2017-07-30 20:45:12
wgdd(月只蓝代发): 金币+20, 感谢帮助！ 2017-08-03 12:14:40
引用回帖:
5楼: Originally posted by wgdd at 2017-07-30 13:35:20
仿照我上面贴出来的论文中的证明过程，可以用留数定理证明出来，当k的虚部为0，解析解就是上面我写的那个的吧。
十分感谢~...

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

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

» 本帖已获得的红花（最新10朵）

wgdd
ส็็็็็็็็็็็็็็็็็็็็
6楼2017-07-30 16:31:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wgdd

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

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

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

如果可以的话，能给些建议吗？应该如何改写，哪儿有类似的程序呢？可以学习一下。
赞一下 回复此楼
7楼2017-07-30 16:35:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cooooldog

铁杆木虫 (著名写手)

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

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

如果能够通用，mathworks 大概率会集成到某个现成的函数上
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

wgdd
ส็็็็็็็็็็็็็็็็็็็็
8楼2017-07-30 17:14:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wgdd

木虫 (正式写手)
送红花一朵
引用回帖:
8楼: Originally posted by cooooldog at 2017-07-30 17:14:59
通用的解决方案比较难，大部分只能针对特定的问题找特定的解法。

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

谢谢，原来这个问题这么复杂~
赞一下 回复此楼
9楼2017-07-30 17:49:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

uvwwawc33
10楼2017-08-17 17:24:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wgdd 的主题更新
361/4返回列表1234下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭