这种类型的二重积分该如何实现?

悠若谷

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

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

rocwoods

这个积分是发散的，可以用放缩法证明。
事实上，原始积分即使放在maple下用符号解，似乎也很难求解出来。一个可行的办法就是用放缩法证明：
分析被平方的积分项，它应该大于

1. 
   
2. int(exp(-1)/(y^2+x^2),x,-1,1)^2
   

复制代码

于是原始积分整体大于

1. 
   
2. 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. 
   
2. syms x y
   
3. a = int(2*y*exp(-y^2)*int(exp(-1)/(y^2+x^2),x,-1,1)^2,0,1)
   
4. double(a)
   
5. Warning: Explicit integral could not be found.
   
6. > In sym.int at 58
   
7. a =
   
8. int(686231412107798092534547949169/633825300114114700748351602688/y*exp(-y^2)*atan(1/y)^2,y = 0 .. 1)
   
9. ans =
   
10.    Inf
    

复制代码

可见结果是无穷的，于是原始积分是无穷的。
所以原始积分不能从0开始积分，下面用我以前给出的方法数值求解下：积分下限分为0.2，0.1，0.01，0.001，0.0001
代码如下：

1. 
   
2. tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.2,1),toc
   
3. f1 =
   
4.    17.9174
   
5. Elapsed time is 0.110964 seconds.
   
6. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.1,1),toc
   
7. f1 =
   
8.    53.9235
   
9. Elapsed time is 0.335609 seconds.
   
10. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.01,1),toc
    
11. f1 =
    
12.   891.1743
    
13. Elapsed time is 1.764475 seconds.
    
14. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.001,1),toc
    
15. f1 =
    
16.   9.7203e+003
    
17. Elapsed time is 8.466365 seconds.
    
18. >> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.0001,1),toc
    
19. f1 =
    
20.   9.8493e+004
    
21. Elapsed time is 21.919113 seconds.
    

复制代码

随着下限逐渐靠近0，所需计算时间也会逐渐增长，这是因为积分是发散的，下限越靠近0，积分步长会分得越细，以期获得足够的精度，当然计算量也会上去。
楼主酌情让y的积分下限取个合适的值吧。

rocwoods

不好意思，上述积分丢掉了2*y*exp(-y^2)项，现在补上：

1. 
   
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.2,1),toc
   
3. 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
   
4. f1 =
   
5.   10.213463739161501
   
6. Elapsed time is 0.056856 seconds.
   
7. f1 =
   
8.   19.868526412761902
   
9. Elapsed time is 0.245491 seconds.
   
10. >> 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
    
11. f1 =
    
12.   61.334993990451785
    
13. Elapsed time is 1.103514 seconds.
    
14. >> 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
    
15. f1 =
    
16.     1.063674743047246e+002
    
17. Elapsed time is 2.785133 seconds.
    
18. >> 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
    
19. f1 =
    
20.     1.517765989598245e+002
    
21. Elapsed time is 5.657963 seconds.
    

复制代码

悠若谷

看到那边论坛的回复了，谢谢您把发散这个问题验证出来了，真的真的很有帮助哪~

悠若谷

还想再问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);