请教数值积分

luobinhan001

我碰到了一个二重奇异振荡积分，做了好久了也没有弄出来，希望各位高手能帮个忙。
由于公式不方便编辑，请各位看附件。
这是一个傅立叶逆变换的一个积分。
方法1:
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=dblquad(@(q1,q2)F,0,infinte,0,infinite)

方法2： 画出被积函数的图象，确定积分限。但是结果不对
alph1=0.5196; alph2=0.9; x=0; y=0.5; z=0;
q1=0:1:500; q2=0:1:500;
for i=1:501,  for j=1:501
        n1(i,j)=sqrt((1-alph1^2)*q1(i)^2+q2(j)^2);
        n2(i,j)=sqrt((1-alph2^2)*q1(i)^2+q2(j)^2);
        B(i,j)=(q1(i)^2+q2(j)^2+n2(i,j)^2)^2-4*n1(i,j)*n2(i,j)*(q1(i)^2+q2(j)^2);
        F(i,j)=n1(i,j)*(q1(i)^2+q2(j)^2+n2(i,j)^2)-2*n1(i,j)*(q1(i)^2+q2(j)^2);
        E(i,j)=cos(q1(i)*x)*cos(q2(j)*y);
        Y(i,j)=F(i,j)*E(i,j)/B(i,j);
end; end
surf(q1,q2,Y); hold on; grid on
xlabel('q1'); ylabel('q2'); zlabel('y'); title('被积函数图像')
f=dblquad(@(q1,q2)F,0,500,0,500)

方法3：采用变步长积分
alph1=0.5196; alph2=0.9; x=0; y=0.5; z=0;
%给数组q1,q2分段赋值，
for i=1:1000                         %步长为0.005
    q1(i)=i*5e-3;                    %q1(1000)=5
    q2(i)=i*5e-3;                    %q2(1000)=5
end
for i=1001:1500                      %步长为0.05
    q1(i)=q2(1000)+(i-1000)*5e-2;    %q1(1500)=5+25=30
    q2(i)=q2(1000)+(i-1000)*5e-2;    %q2(1500)=5+25=30
end
for i=1501:2000                      %步长为0.5
    q1(i)=q2(1500)+(i-1500)*5e-1;    %q1(2000)=30+250=280
    q2(i)=q2(1500)+(i-1500)*5e-1;    %q2(2000)=30+250=280
end
for i=2001:2100                      %步长为5
    q1(i)=q2(2000)+(i-2000)*5;       %q2(2000)=280+500=780
    q2(i)=q2(2000)+(i-2000)*5;       %q2(2000)=280+500=780
end
   
for i=1:2100,   for j=1:2100
        n1(i,j)=sqrt((1-alph1^2)*q1(i)^2+q2(j)^2);
        n2(i,j)=sqrt((1-alph2^2)*q1(i)^2+q2(j)^2);
        B(i,j)=(q1(i)^2+q2(j)^2+n2(i,j)^2)^2-4*n1(i,j)*n2(i,j)*(q1(i)^2+q2(j)^2);
        F(i,j)=n1(i,j)*(q1(i)^2+q2(j)^2+n2(i,j)^2)-2*n1(i,j)*(q1(i)^2+q2(j)^2);
        E(i,j)=cos(q1(i)*x)*cos(q2(j)*y);
        Y(i,j)=F(i,j)*E(i,j)/B(i,j);
end; end
surf(q1,q2,Y); hold on; grid on
xlabel('q1'); ylabel('q2'); zlabel('y'); title('被积函数图像')
f=dblquad(@(q1,q2)Y,1e-3,556,1e-3,556);
final=f/(pi^2*2e7)

发现结果数量极都不对。
想法：被积函数在点（0，0）处，分母是0的4次方，分子是0的3次方，应该是无穷大，类似于对1/x 积分的感觉，我觉得解析解应该是无穷，也就是对于我这个积分，如果子步越小，得出来的值应该越大。要想给出一个满意的数值解，我觉得不大可能。
希望各位高手给于解答。

[ 本帖最后由 ChaChing 于 2010-8-3 20:45 编辑 ]

xjzuo

画了一下被积函数的图形,发现有点奇异.
所以建议改用离散求和的形式求积分.

luobinhan001

这样做对不对？
alph1=0.5196;
alph2=0.9;
x=0;
y=0.5;
z=0;
q1=-500:10:500;
q2=-500:10:500;
for i=1:101
    for j=1:101
        n1(i,j)=sqrt((1-alph1^2)*q1(i)^2+q2(j)^2);
        n2(i,j)=sqrt((1-alph2^2)*q1(i)^2+q2(j)^2);
        B(i,j)=(q1(i)^2+q2(j)^2+n2(i,j)^2)^2-4*n1(i,j)*n2(i,j)*(q1(i)^2+q2(j)^2);
        F(i,j)=n1(i,j)*(q1(i)^2+q2(j)^2+n2(i,j)^2)-2*n1(i,j)*(q1(i)^2+q2(j)^2);
        E(i,j)=cos(q1(i)*x)*cos(q2(j)*y);
        Y(i,j)=F(i,j)*E(i,j)/B(i,j);
    end
end
surf(q1,q2,Y)
hold on
grid on
xlabel('q1')
ylabel('q2')
zlabel('y')
title('被积函数图像')


然后在matlab的窗口
dblquad(@(q1,q2)Y,-500,500,-500,500)
ans=393.6004

其实从图中可以看出q1在正负300左右非常小，接近于零。如果改变积分限
dblquad(@(q1,q2)Y,-300,300,-500,500)
ans=236.1131
差距怎么会这么大呢？

另外如果我把被积函数的变量区间取得更大一点的话，（-1000，1000）
理论上结果应该相差不大。但是matlab给出的结果让人无法接受。
alph1=0.5196;
alph2=0.9;
x=0;
y=0.5;
z=0;
q1=-1000:10:1000;
q2=-1000:10:1000;
for i=1:201
    for j=1:201
        n1(i,j)=sqrt((1-alph1^2)*q1(i)^2+q2(j)^2);
        n2(i,j)=sqrt((1-alph2^2)*q1(i)^2+q2(j)^2);
        B(i,j)=(q1(i)^2+q2(j)^2+n2(i,j)^2)^2-4*n1(i,j)*n2(i,j)*(q1(i)^2+q2(j)^2);
        F(i,j)=n1(i,j)*(q1(i)^2+q2(j)^2+n2(i,j)^2)-2*n1(i,j)*(q1(i)^2+q2(j)^2);
        E(i,j)=cos(q1(i)*x)*cos(q2(j)*y);
        Y(i,j)=F(i,j)*E(i,j)/B(i,j);
    end
end
surf(q1,q2,Y)
hold on
grid on
xlabel('q1')
ylabel('q2')
zlabel('y')
title('被积函数图像')


在matlab窗口
dblquad(@(q1,q2)Y,-1000,1000,-1000,1000)
ans=-2.8896*e3

不仅数量级都变了，而且符号都相反了。

是不是对于这种奇异积分不能这样做？

luobinhan001

求一个类似于 （x+y）/(xy)这样函数的积分，积分区间为（0，1，0，1）
在论文中遇到一个无穷限的广义积分。
syms x y
f=(x+y)/(x*y);
dblquad(@(x,y)f,0,1,0,1)


错误提示为：
??? Error using ===> innerintegral
Inputs to innerintegral must be floats, namely single or double.
Error in ==> dblquad>innerintegral at 79
Q = zeros(size(y), superiorfloat(fcl, xmax, y));
Error in ==> quad at 63
y = f(x, varargin{:});
Error in ==> dblquad at 58
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...
Error in ==> shiyan at 3
dblquad(@(x,y)f,0,1,0,1)



希望各位高手发表意见，不胜感激。

[ 本帖最后由 ChaChing 于 2010-8-3 20:47 编辑 ]

xjzuo

将你的原问题贴一下.
你举的这个例子本身就发散,并不是广义无穷积分.

luobinhan001

二重奇异积分求助

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

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的用法不熟悉吧

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

xjzuo

这个问题好象已经讨论过------奇异积分用离散求和较好,dblquad一般是失效的.

luobinhan001

谢谢楼上。现将问题详细描述如下：

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)
这要如何处理呢？

导师催得紧，都快要疯了。


[ 本帖最后由 ChaChing 于 2010-8-3 20:33 编辑 ]

xjzuo

晕!你竟然用好几个马甲在这里重复问一个问题!----真不知道这样是否有违版规?

另:恐怖的是,你竟然还在仿真论坛上自称"小女子"!

[ 本帖最后由 xjzuo 于 2007-5-17 23:31 编辑 ]

luobinhan001

这是我同学的一个问题，她找我帮忙，我也搞不定。所以就发上来了，希望人多力量大。如有违规之处还请见谅。对不起