声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1601|回复: 4

[综合讨论] 这种类型的二重积分该如何实现?

[复制链接]
发表于 2009-9-14 13:05 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
积分.jpg
积分方程如图所示
用哪个函数实现比较好呢?
dblquad?
恳请各位前辈帮忙看看!

[ 本帖最后由 悠若谷 于 2009-9-14 13:06 编辑 ]

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2009-9-16 01:28 | 显示全部楼层
这个积分是发散的,可以用放缩法证明。
事实上,原始积分即使放在maple下用符号解,似乎也很难求解出来。一个可行的办法就是用放缩法证明:
分析被平方的积分项,它应该大于

  1. int(exp(-1)/(y^2+x^2),x,-1,1)^2
复制代码
于是原始积分整体大于

  1. int(2*y*exp(-y^2)*int(exp(-1)/(y^2+x^2),x,-1,1)^2,0,1)
复制代码
在MATLAB2008a下运行如下代码:(之所以选择2008a,是因为这是MATLAB和Maple结合的最后一个版本,以后的版本符号计算内核换了mupad,下列代码是得不出结果的,可见mupad比maple还是弱多了)

  1. syms x y
  2. a = int(2*y*exp(-y^2)*int(exp(-1)/(y^2+x^2),x,-1,1)^2,0,1)
  3. double(a)
  4. Warning: Explicit integral could not be found.
  5. > In sym.int at 58
  6. a =
  7. int(686231412107798092534547949169/633825300114114700748351602688/y*exp(-y^2)*atan(1/y)^2,y = 0 .. 1)
  8. ans =
  9.    Inf
复制代码
可见结果是无穷的,于是原始积分是无穷的。
所以原始积分不能从0开始积分,下面用我以前给出的方法数值求解下:积分下限分为0.2,0.1,0.01,0.001,0.0001
代码如下:

  1. tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.2,1),toc
  2. f1 =
  3.    17.9174
  4. Elapsed time is 0.110964 seconds.
  5. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.1,1),toc
  6. f1 =
  7.    53.9235
  8. Elapsed time is 0.335609 seconds.
  9. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.01,1),toc
  10. f1 =
  11.   891.1743
  12. Elapsed time is 1.764475 seconds.
  13. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.001,1),toc
  14. f1 =
  15.   9.7203e+003
  16. Elapsed time is 8.466365 seconds.
  17. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.0001,1),toc
  18. f1 =
  19.   9.8493e+004
  20. Elapsed time is 21.919113 seconds.
复制代码
随着下限逐渐靠近0,所需计算时间也会逐渐增长,这是因为积分是发散的,下限越靠近0,积分步长会分得越细,以期获得足够的精度,当然计算量也会上去。
楼主酌情让y的积分下限取个合适的值吧。

评分

1

查看全部评分

发表于 2009-9-17 22:11 | 显示全部楼层
不好意思,上述积分丢掉了2*y*exp(-y^2)项,现在补上:

  1. tic;f1 = quadl(@(y) 2*y.*exp(-y.^2).*arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.2,1),toc
  2. tic;f1 = quadl(@(y) 2*y.*exp(-y.^2).*arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.1,1),toc
  3. f1 =
  4.   10.213463739161501
  5. Elapsed time is 0.056856 seconds.
  6. f1 =
  7.   19.868526412761902
  8. Elapsed time is 0.245491 seconds.
  9. >> tic;f1 = quadl(@(y) 2*y.*exp(-y.^2).*arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.01,1),toc
  10. f1 =
  11.   61.334993990451785
  12. Elapsed time is 1.103514 seconds.
  13. >> tic;f1 = quadl(@(y) 2*y.*exp(-y.^2).*arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.001,1),toc
  14. f1 =
  15.     1.063674743047246e+002
  16. Elapsed time is 2.785133 seconds.
  17. >> tic;f1 = quadl(@(y) 2*y.*exp(-y.^2).*arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.0001,1),toc
  18. f1 =
  19.     1.517765989598245e+002
  20. Elapsed time is 5.657963 seconds.
复制代码

评分

1

查看全部评分

 楼主| 发表于 2009-9-18 20:49 | 显示全部楼层

回复 板凳 rocwoods 的帖子

看到那边论坛的回复了,谢谢您把发散这个问题验证出来了,真的真的很有帮助哪~
 楼主| 发表于 2009-9-18 22:52 | 显示全部楼层

回复 板凳 rocwoods 的帖子

还想再问rocwoods前辈一个问题,
我把这种类型的积分跟fzero函数结合来解方程的时候
出现??? Error using ==> fzero at 317
FZERO cannot continue because user supplied function_handle ==>
@(N)quadl(@(y)k*sqrt((ns+dN*y/(pi*D).*exp(-(y/D).^2).*arrayfun(@(s)quadl(@(s)exp(-s.^2)./(s.^2+(y/D).^2),-w/D,w/D),s)).^2-N^2),1e-10,3e-6)-pi/4-atan(sqrt((N.^2-nc^2)/(N1^2-N.^2))+(k^2*N1*dN*(-2/pi)*(exp(-(w/D)^2)/w+pi^0.5/D*erf(w/D)))/(2*(k^2*N1^2-k^2*N.^2)^1.5))
failed with the error below.
不知道该怎么处理
敲入的代码如下,方程非常复杂
ns=1.520167;
nf=1.521484;
k=2*pi/0.6328e-6;
nc=1;
w=2.5e-6;
D=3e-6;
N1=fzero(@(y) 2*w*k*sqrt(nf^2-y^2)-pi+2*atan(ns^2/(nf^2)*sqrt((nf^2-y^2)/(y^2-ns^2))),[1.520168,1.521483]);
dN=N1-ns;
N=fzero(@(N) quadl(@(y) k*sqrt((ns+dN*y/(pi*D).*exp(-(y/D).^2).*arrayfun(@(s) quadl(@(s) exp(-s.^2)./(s.^2+(y/D).^2),-w/D,w/D),s)).^2-N^2),1e-10,3e-6)-pi/4-atan(sqrt((N.^2-nc^2)/(N1^2-N.^2))+(k^2*N1*dN*(-2/pi)*(exp(-(w/D)^2)/w+pi^0.5/D*erf(w/D)))/(2*(k^2*N1^2-k^2*N.^2)^1.5)),1.52);
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 19:57 , Processed in 0.077965 second(s), 29 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表