这种类型的二重积分该如何实现?
积分方程如图所示
用哪个函数实现比较好呢？
dblquad？
恳请各位前辈帮忙看看！
[ 本帖最后由 悠若谷 于 2009-9-14 13:06 编辑 ] 这个积分是发散的，可以用放缩法证明。
事实上，原始积分即使放在maple下用符号解，似乎也很难求解出来。一个可行的办法就是用放缩法证明：
分析被平方的积分项，它应该大于
int(exp(-1)/(y^2+x^2),x,-1,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还是弱多了）
syms x y
a = int(2*y*exp(-y^2)*int(exp(-1)/(y^2+x^2),x,-1,1)^2,0,1)
double(a)
Warning: Explicit integral could not be found.
> In sym.int at 58
a =
int(686231412107798092534547949169/633825300114114700748351602688/y*exp(-y^2)*atan(1/y)^2,y = 0 .. 1)
ans =
Inf
可见结果是无穷的，于是原始积分是无穷的。
所以原始积分不能从0开始积分，下面用我以前给出的方法数值求解下：积分下限分为0.2，0.1，0.01，0.001，0.0001
代码如下：
tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.2,1),toc
f1 =
17.9174
Elapsed time is 0.110964 seconds.
>> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.1,1),toc
f1 =
53.9235
Elapsed time is 0.335609 seconds.
>> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.01,1),toc
f1 =
891.1743
Elapsed time is 1.764475 seconds.
>> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.001,1),toc
f1 =
9.7203e+003
Elapsed time is 8.466365 seconds.
>> tic;f1 = quadl(@(y) arrayfun(@(y)quadl(@(x) exp(-x.^2)./(y.^2+x.^2),-1,1),y).^2,0.0001,1),toc
f1 =
9.8493e+004
Elapsed time is 21.919113 seconds.
随着下限逐渐靠近0，所需计算时间也会逐渐增长，这是因为积分是发散的，下限越靠近0，积分步长会分得越细，以期获得足够的精度，当然计算量也会上去。
楼主酌情让y的积分下限取个合适的值吧。 不好意思，上述积分丢掉了2*y*exp(-y^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
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
f1 =
10.213463739161501
Elapsed time is 0.056856 seconds.
f1 =
19.868526412761902
Elapsed time is 0.245491 seconds.
>> 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
f1 =
61.334993990451785
Elapsed time is 1.103514 seconds.
>> 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
f1 =
1.063674743047246e+002
Elapsed time is 2.785133 seconds.
>> 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
f1 =
1.517765989598245e+002
Elapsed time is 5.657963 seconds.
回复 板凳 rocwoods 的帖子
看到那边论坛的回复了，谢谢您把发散这个问题验证出来了，真的真的很有帮助哪~回复 板凳 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))),);
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);
页:
[1]