二重奇异积分求助

shangguanwaner

小女子在做毕业论文的时候碰上了一个二重积分，做了好久都求不出来，求各位高手帮忙。不甚感激。

syms q1 q2
alph1=0.5196;
alph2=0.9;
x=0;
y=0.5;
z=0;
n1=sqrt((1-alph1^2)*q1^2+q2^2);
n2=sqrt((1-alph2^2)*q1^2+q2^2);
B=(q1^2+q2^2+n2^2)^2-4*n1*n2*(q1^2+q2^2);
F=(n1*(q1^2+q2^2+n2^2)*exp(-n1*z)-2*n1*(q1^2+q2^2)*exp(-n2*z))*cos(q1*x)*cos(q2*y)/B;



F是被积函数 ，积分区间为（-无穷，＋无穷），（-无穷，＋无穷）
当然积分区间也可以不必为无穷，q1,q2在1000左右的时候，F的值已经很小了。
参考了http://www.simwe.com/forum/archiver/tid-490651.html后。
1我用了“午夜流星”的方法
g = dblquad(inline(F), 0+eps, 1000, 0+eps, 1000);
D=g/(4*pi^2*20)
出现错误：好像是提示分母B＝0 。我不解的是为什么“tooth52”的函数也有分母为零的时候，为什么她的能算出结果来呢？而我的却不行。


2我把zzgrnr的命令流复制到matlab里面运行失败，提示Undefined function or variable "wfxy".
由于没有看懂程序也不敢随便改里面的WFXY。


最后也试了一下bainhome提供的工具箱NIT，M文件总是出错，可那是对inline的用法不熟悉吧

请各位大侠无论如何帮个忙。我都心力憔悴了。
在线等。

rocwoods

请参考下面的例子。
http://bbs.matwav.com/post/view? ... 4&sty=1&tpg=2&age=0

shangguanwaner

看了楼上的连接，感觉和我这个问题没有什么相似的地方。麻烦楼上GG说得详细一点好吗

rocwoods

直接用匿名函数@表达式比较麻烦。干脆用下面的形式吧：nested function
function g=shangguanwaner(alph1,alph2,x,y,z)
alph1=0.5196;
alph2=0.9;
x=0;
y=0.5;
z=0;
function F=jifen(q1,q2)
n1=sqrt((1-alph1^2).*q1.^2+q2.^2);
n2=sqrt((1-alph2^2).*q1.^2+q2.^2);
B=(q1.^2+q2.^2+n2.^2).^2-4*n1.*n2.*(q1.^2+q2.^2);
F=(n1.*(q1.^2+q2.^2+n2.^2).*exp(-n1.*z)-2.*n1.*(q1.^2+q2.^2).*exp(-n2.*z)).*cos(q1.*x).*cos(q2.*y)./(B+eps*(B==0));
end
g=dblquad(@jifen,-600,600,-600,600);
end

正确把积分表达式写出来应该没问题。注意要写成如上面点乘的形式。最后(B+eps*(B==0))是为了防止除以0。
把上面文件保存后运行g=shangguanwaner(0.5196,0.90,0,0.5,0)即可

shangguanwaner

谢谢rocwoods。我运行你的程序的以后，得到的结果当y=0.5的时候，结果： 1.119499034512188e+002
是比较可靠的。但是我令y=1.0时，得到的结果：-9.603938278187914e+002。这两个结果感觉不应该是异号的。而且应该是随着y的增大，结果越来越小。每次运行时候总是弹出：Warning: Minimum step size reached; singularity possible. 另外想问个问题：出现这个警告的时候怎么评估结果的可靠性呢？这个被积函数确实有奇异性。当q1=q2=0时，分母是0的4次方，而分子是0的3次方。
另外你的这种表示方法确实感到很方便(B+eps*(B==0))
再一次谢谢rocwoods

rocwoods

评估结果的可靠性，就需要详细了解dblquad所采用的算法了。我现在是没有精力研究这个算法了。此外，没有手写的公式看着也太费劲了，你可以贴个公式图上来，大家帮你分析分析。
建议，不知道成不成，可以避开0点积分。将q1，q2的积分区域都分成(-600,-realmin)与(realmin,600)的和。
明天要出差了，祝你好运。

shangguanwaner

太感谢rocwoods了.刚才用了你的方法，得到的结果要好很多
function g=shangguanwaner(alph1,alph2,x,y,z)
alph1=0.5196;
alph2=0.9;
x=0;
y=0.5;
z=0;
function F=jifen(q1,q2)
n1=sqrt((1-alph1^2).*q1.^2+q2.^2);
n2=sqrt((1-alph2^2).*q1.^2+q2.^2);
B=(q1.^2+q2.^2+n2.^2).^2-4*n1.*n2.*(q1.^2+q2.^2);
F=(n1.*(q1.^2+q2.^2+n2.^2).*exp(-n1.*z)-2.*n1.*(q1.^2+q2.^2).*exp(-n2.*z)).*cos(q1.*x).*cos(q2.*y)./B;
end
g=dblquad(@jifen,realmin,600,realmin,600);
G=g/(pi^2*20)
end



当y=0.5,1.0,1.5时，G＝0.1418，0.0750，0.0525    这三个结果与其他文献结果都非常接近

然而当y=2.0，2.5时,G=-0.1296，-0.0176  这两个结果就不对了
理论上是随着y的增大，G越来越小接近于0。
这是不是函数的振荡性造成的？被积函数含有cos(q1.*x).cos(q2.*y)
这要如何处理呢？

:'( 不会传公式的图片。上传附件说是图片太大了，格式改成gif,jpg,jpeg都没有减小文件:'(

:D 好不容易上传成功了

shangguanwaner

晕倒，图片怎么变成这样了。
再试一次