嵌套二重积分的插值解法和RBF网络解法

rocwoods

本贴探讨嵌套二重积分的一种插值解法。希望抛砖引玉，大家踊跃讨论。以http://www.simwe.com/forum/thread-789115-1-1.html为例。
      上述帖子中的二重积分内部还有两个二重积分的平方。按通常办法解起来，计算量相当大，以我在那个帖子里给出的均匀划分网格法为例，虽然该方法将积分区域划分成均匀的n*n个网格，以每个网格中心的函数值乘以网格的面积作为函数在这个网格上积分的近似，由于内部函数还要计算两个二重积分的平方，因此这一方法要计算2*n^2个二重积分的平方。计算量还是太大。
      本贴的想法是将插值法和均匀划分网格法结合起来，为了求一个复杂的连续性、光滑性比较好的函数（譬如上面帖子链接里的函数）在一个区域上的积分，首先粗划分积分区域（生成比较粗的网格），计算函数在粗网格节点上的值，然后在上面网格基础上细划分网格，根据粗网格节点上的函数值利用interp2或者griddata函数插值生成细网格节点上的函数值。然后再乘以细网格面积并求和得到整个函数在积分区域积分的近似值。以下是以上面提到的积分为例子的代码：

1. function [INT,TfI]=alber2008i(inf_u,inf_v,mm,nn)
   
2. %inf_u,u无穷上限的近似
   
3. %inf_v，v无穷上限的近似,
   
4. %mm,u,v积分限细等分个数(目的是插值生成足够多的Tf_Integrand函数点)
   
5. %nn,u,v积分限粗等分个数(目的是生成插值的骨架点)
   
6. [UI,VI]=meshgrid( (inf_u/mm:inf_u/mm:inf_u),(inf_v/mm:inf_v/mm:inf_v)) ;%细网格数据
   
7. [U,V]=meshgrid( linspace(inf_u/mm,inf_u,nn),linspace(inf_v/mm,inf_v,nn)) ;%粗网格数据
   
8. Tf=zeros(nn);
   
9. for k=1:nn^2   
   
10.         Tf(k)=Tf_Integrand( U(k), V(k) ); %生成粗网格节点函数值
    
11. end
    
12. %插值并求近似积分
    
13. TfI = interp2(U,V,Tf,UI,VI);
    
14. INT=TfI*inf_u*inf_v/(mm^2);
    
15. INT=sum(INT(:));
    
16. function z = Tf_Integrand(u,v)
    
17. m=1;
    
18. n=1;
    
19. a = 1;
    
20. b = 1;
    
21. km = 2*m*pi/a;
    
22. kn = 2*n*pi/b;
    
23. f_rel = @(x,y) (1-cos(km*x)).*(1-cos(kn*y)).*cos(u.*x+v.*y);
    
24. f_img = @(x,y) (1-cos(km*x)).*(1-cos(kn*y)).*sin(u.*x+v.*y);
    
25. tmp_rel = dblquad(f_rel, 0,a, 0,b);
    
26. tmp_img = dblquad(f_img, 0,a, 0,b);
    
27. z = (tmp_rel.^2+tmp_img.^2)./sqrt(u.^2+v.^2);
    
28. end
    
29. end

复制代码

1. 划分网格法：分成25*25个网格，计算1250个二重积分的平方：所用时间结果：
   
2. tic,INT=alber2008(40,40,25);toc
   
3. Elapsed time is 188.849775 seconds.
   
4. INT =
   
5.     6.8267
   
6. 插值法：粗分成15*15个网格，计算450个二重积分的平方，然后根据这225个函数点插值生成50*50个函数值，最后计算近似积分，所用时间和结果：
   
7. tic;[INTi,TfI]=alber2008i(40,40,50,15);toc
   
8. Elapsed time is 68.853558 seconds.
   
9. INTi =
   
10.     6.0241

复制代码

[ 本帖最后由 rocwoods 于 2007-7-28 23:54 编辑 ]

taohe

挺好的想法。能有效的提高求解速度。这个是否类似于梯度法求解定积分？

从这个例子，我学到解决实际问题的一个很好的办法就是像这样，在掌握基本原则的基础上，针对问题本身，选择有效的方法。谢谢。

可惜的是，原帖发问的同学并没有跟上讨论。。。(我突然想起来了有个帖子似乎是发难说这里为什么没有人回答问题什么什么的，有点奇怪）

rocwoods

我不太了解梯度法求解定积分具体实现过程。不管那个同学讨论不讨论，还得感谢他提出的这个问题。今天又想了想这个问题，发现插值方法虽然比最初的方法快，可是精度是个问题，如何同时保证精度和速度是我们努力的方向。
       今天突然想起以前做过的RBF网络逼近函数的一个例子。RBF网络是人工神经网络的一种，有”局部感知场“美誉，可以任意逼近光滑性比较好的函数。下面有我之前关于它的一个讨论。
http://bbs.matwav.com/post/view?bid=12&id=560285&tpg=1&ppg=1&sty=1&age=0#560285
除了上面帖子链接里给出的RBF网络对正弦函数这个一维函数逼近的例子，下面再给一个RBF网络对MATLAB里peaks这个二维函数逼近的例子。大家可以看看它比插值法的逼近精度高多少!

1. [X,Y] = meshgrid(-3:.25:3);
   
2. Z = peaks(X,Y);%生成网络训练样本
   
3. input=[X(:)';Y(:)'];%网络输入
   
4. output=[Z(:)'];%网络输出
   
5. net=newrb(input,output,0.00001,0.5);%生成peaks的逼近网络
   
6. [XI,YI] = meshgrid(-3:.125:3);
   
7. Znet=sim(net,[XI(:)';YI(:)']);%网络给出的peaks函数值
   
8. Znet=reshape(Znet,49,49);
   
9. 
   
10. Zpeak=peaks(XI,YI);%真实的peaks函数值
    
11. err1=abs(Znet-Zpeak);%网络逼近值和真实值的绝对误差
    
12. Zinter = interp2(X,Y,Z,XI,YI);%插值给出的peaks函数值
    
13. err2=abs(Zinter-Zpeak);%插值逼近值和真实值的绝对误差
    
14. sum(err1(:))
    
15. sum(err2(:))

复制代码

从上面可以看出RBF网络的逼近能力远高于插值法的逼近能力。于是就有了下面的方法（目前我能想到的最好的方法）：

1. function [INT,TfInet]=rbf_calc(inf_u,inf_v,nn)
   
2. %inf_u,u无穷上限的近似
   
3. %inf_v，v无穷上限的近似,
   
4. %nn，u,v积分限区间的等分个数
   
5. 
   
6. [UI,VI]=meshgrid(linspace(0,inf_u,nn+1),linspace(0,inf_v,nn+1));
   
7. UI=(UI(1:end-1,1:end-1)+UI(2:end,2:end))/2;
   
8. VI=(VI(1:end-1,1:end-1)+VI(2:end,2:end))/2;%生成nn*nn网格的中心点
   
9. 
   
10. [U,V]=meshgrid(linspace(UI(1,1),UI(1,end),nn/1.5),linspace(VI(1,1),VI(end,1),nn/1.5));%生成网络训练样本网格,样本网格均匀分布在nn*nn网格区域内，密度为UI,VI
    
11. %网格的1/1.5
    
12. Tf=zeros(size(U));
    
13. for k=1:length(Tf)^2   
    
14.         Tf(k)=Tf_Integrand( U(k),V(k));   
    
15. end
    
16. %训练网络并给出足够多的Tf_Integrand函数点并计算积分
    
17. input=[U(:)';V(:)'];%网络输入
    
18. output=[Tf(:)'];%网络输出
    
19. net=newrb(input,output,0.00001,3*inf_u/nn);%生成逼近网络
    
20. 
    
21. 
    
22. TfInet=sim(net,[UI(:)';VI(:)']);%网络给出的函数值
    
23. TfInet=reshape(TfInet,nn,nn);
    
24. 
    
25. INT=TfInet*inf_u*inf_v/(nn^2);
    
26. INT=sum(INT(:));
    
27. 
    
28. function z = Tf_Integrand(u,v)
    
29. m=1;
    
30. n=1;
    
31. a = 1;
    
32. b = 1;
    
33. km = 2*m*pi/a;
    
34. kn = 2*n*pi/b;
    
35. 
    
36. f_rel = @(x,y) (1-cos(km*x)).*(1-cos(kn*y)).*cos(u.*x+v.*y);
    
37. f_img = @(x,y) (1-cos(km*x)).*(1-cos(kn*y)).*sin(u.*x+v.*y);
    
38. tmp_rel = dblquad(f_rel, 0,a, 0,b);
    
39. tmp_img = dblquad(f_img, 0,a, 0,b);
    
40. 
    
41. z = (tmp_rel.^2+tmp_img.^2)./sqrt(u.^2+v.^2);
    
42. end
    
43. end

复制代码

以下是和第一种方法的对比：

1. tic,[INT,intt]=alber2008(40,40,30);toc
   
2. Elapsed time is 260.640272 seconds.
   
3. >> INT
   
4. INT =
   
5.     6.9354
   
6. 
   
7. tic,[INTnet,TfInet]=rbf_calc(40,40,30);toc
   
8. NEWRB, neurons = 0, SSE = 0.45778
   
9. NEWRB, neurons = 25, SSE = 0.047086
   
10. NEWRB, neurons = 50, SSE = 0.00730593
    
11. NEWRB, neurons = 75, SSE = 0.00114306
    
12. NEWRB, neurons = 100, SSE = 0.000899809
    
13. NEWRB, neurons = 125, SSE = 0.000185916
    
14. NEWRB, neurons = 150, SSE = 3.66627e-005
    
15. Elapsed time is 121.612190 seconds.
    
16. >> INTnet
    
17. INTnet =
    
18.     7.0517

复制代码

可见从精度和速度上都兼顾了!事实上，训练好网络后，完全可以再更加多的细分网格，因为net是Tf_Integrand函数的近似，给定一组特定的(u,v),用sim(net,[u;v])可以非常迅速求出Tf_Integrand函数的近似值。这样求最后的结果就不用担心网格太密的时候计算量大了。有兴趣的朋友可以试试。

[ 本帖最后由 rocwoods 于 2007-7-28 23:52 编辑 ]

bainhome

请教u和v这两个变量两位是如何处理的？taohe的程序我加断点调试发现利用dblquad计算中，u和v均等于0.1，也就是调用中间函数时，嵌套的那两个dblquad所使用的u和v不是未知变量，而是直接等于如下代码中的第二个和第四个参数0.1、0.1

1. tic; f2=dblquad(@Tf_Integrand, 0.1, 1, 0.1, 1), toc

复制代码

请教下调用的原因为何？
另外，此类矩形区域二重积分，似乎也可以通过变量代换到[-1,1]范围内利用Gauss积分进行计算，虽然Gauss积分法只适合于多项式被积函数，不过参考数值分析书中所述，似可通过增加积分点的个数得到改善。目前的困难是变量替换我还没太大把握,计算结果差别很大。给出我的代码，其原理多见于各种有限元和数值分析书中不再赘述：

1. % ****************** 两 * 点 * Gauss 积 * 分 * 计 * 算 * 程 * 序 ********************
   
2. clc
   
3. m=1;
   
4. n=1;
   
5. a = 1;
   
6. b = 1;
   
7. km = 2*m*pi/a;
   
8. kn = 2*n*pi/b;
   
9. %＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
   
10. u=.1;这两句定义不准确，随便填的。
    
11. v=1;
    
12. %＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
    
13. ss=[.5773 .5773 -.5773 -.5773];
    
14. tt=[.5773 -.5773 .5773 -.5773];
    
15. f_rel = @(s,t)(1-cos(km*.5*(s+1)))*(1-cos(kn*.5*(t+1))).*cos(u.*.5*(s+1)+v.*.5*(t+1));
    
16. f_img = @(s,t)(1-cos(km*.5*(s+1)))*(1-cos(kn*.5*(t+1))).*sin(u.*.5*(s+1)+v.*.5*(t+1));
    
17. F1=[];
    
18. F2=[];
    
19. for i=1:4
    
20.     F1=[F1,f_rel(ss(i),tt(i))];
    
21.     F2=[F2,f_img(ss(i),tt(i))];
    
22. end
    
23. F1=sum(F1),F2=sum(F2),
    
24. z = (F1^2+F2^2)/sqrt(u^2+v^2)

复制代码

假如此方法经调整改进可以实现，我有把握将计算时间缩短至rocwoods兄的1/10，but now...:lol

taohe

好方法真不少，都值得我们学习。谢谢分享。

to bainhome：
u，v是被积函数的自变量，如果直接调用dblquad，我们只需给出积分区间：比如u in [0.1 1], v in [0.1 1]。

rocwoods

其实说白了，Tf_Integrand函数是关于u,v的二元函数，只不过这个函数含有两个二重积分的平方，很复杂。因此，三种方法都是基于积分的基本定义来计算的，即细分积分区域，每个小积分区域的面积乘以该区域上的函数值(面积足够小的话可以认为该区域函数值一致)然后再求和。我的第一种方法的思想就是，网格化积分区域后直接计算Tf_Integrand函数在各个网格中心的函数值作为该网格区域的函数值，然后乘以各个网格面积并求和作为积分值得近似。因此这种方法随着积分区域的扩大，保证精度的话计算量会非常大。
第二种方法是用插值的方法去近似Tf_Integrand，只需要计算出Tf_Integrand在一些格点上的值，剩下的格点上的函数值插值得出。这种方法精度很有限。
第三种方法是用RBF神经网络去近似Tf_Integrand，也只需给出Tf_Integrand在一些格点上的值，剩下的格点上的函数值由近似网络给出。RBF神经网络逼近函数的能力很高，因此，我认为该方法优于第二种方法。
     各种数值积分的基本思想都是用简单函数近似代替被积函数，我对简单的理解就是能够非常方便快速的求出函数值。原题中Tf_Integrand这个函数很复杂，之所以复杂是因为计算量很大，为了得到(u,v)上的函数值，要付出计算两个二重积分的平方的代价！这比那些由初等函数复合而成的函数复杂的多。gauss积分法还是基于插值公式的，不同的是，以一重定积分为例，它并不简简单单用插值函数值乘以区间长度并求和作为积分的近似值。因此精度比用插值函数值乘以区间长度并求和要高的多。
      但是当区间长度很长时，直接用Gauss积分法的话，精度也是个问题，以一维情况为例，int(f(x),x,a,b)用Gauss积分法取5个节点时的余式R5=(1/1237732650)*((b-a)/2)^11*f^(10)(kesi)。可想而知，(b-a)/2=20的时候,总误差还是蛮大的。相应于一维情况的5个节点，二维情况就得取5*5个节点。为了获取较高的精度，势必多取节点，在二维情况下节点数就得平方。而每多取一个节点，为了获得节点上的值，就得计算两个二重积分的平方。
      我现在的改进想法就是基于第三种方法。既然已经用RBF网络逼近原函数了，剩下的就是如何利用逼近好的网络采取有效的积分方法(不是简单的用面积乘以函数值)来计算积分。可以用dblquad调用网络来实现，有空试试。

[ 本帖最后由 rocwoods 于 2007-7-30 12:55 编辑 ]

rocwoods

第四种方法：用dblquad调用网络。
为了简单说明这个方法先给一个例子，peaks函数的例子

1. dblquad(@(x,y) peaks(x,y),-3,3,-3,3)
   
2. ans =
   
3.   13.057098664914852

复制代码

dblquad调用网络得出的值：

1. [X,Y] = meshgrid(-3:.25:3);
   
2. Z = peaks(X,Y);%生成网络训练样本
   
3. input=[X(:)';Y(:)'];%网络输入
   
4. output=[Z(:)'];%网络输出
   
5. net=newrb(input,output,0.00001,0.5);%生成peaks的逼近网络
   
6. dblquad(@(x,y) sim(net,[x(:)';y*ones(size(x(:)') )]),-3,3,-3,3)
   
7. NEWRB, neurons = 0, SSE = 1189.19
   
8. NEWRB, neurons = 25, SSE = 7.94428
   
9. NEWRB, neurons = 50, SSE = 0.531099
   
10. NEWRB, neurons = 75, SSE = 0.0676809
    
11. NEWRB, neurons = 100, SSE = 0.00777982
    
12. NEWRB, neurons = 125, SSE = 0.00143653
    
13. NEWRB, neurons = 150, SSE = 0.000216292
    
14. NEWRB, neurons = 175, SSE = 6.17855e-005
    
15. NEWRB, neurons = 200, SSE = 2.02683e-005
    
16. ans =
    
17.   13.057095719412455

复制代码

可以看出精度非常好了。
回到原题：
原题第四种方法代码如下：

1. function INT=rbf_dblquad(inf_u,inf_v,nn)
   
2. %inf_u,u无穷上限的近似
   
3. %inf_v，v无穷上限的近似,
   
4. %nn，u,v积分限区间的等分个数
   
5. 
   
6. [UI,VI]=meshgrid(linspace(0,inf_u,nn+1),linspace(0,inf_v,nn+1));
   
7. UI=(UI(1:end-1,1:end-1)+UI(2:end,2:end))/2;
   
8. VI=(VI(1:end-1,1:end-1)+VI(2:end,2:end))/2;%生成nn*nn网格的中心点
   
9. 
   
10. [U,V]=meshgrid(linspace(UI(1,1),UI(1,end),nn/1.5),linspace(VI(1,1),VI(end,1),nn/1.5));%生成网络训练样本网格,样本网格均匀分布在nn*nn网格区域内，密度为UI,VI
    
11. %网格的1/1.5
    
12. Tf=zeros(size(U));
    
13. for k=1:length(Tf)^2   
    
14.         Tf(k)=Tf_Integrand( U(k),V(k));   
    
15. end
    
16. %训练网络并给出足够多的Tf_Integrand函数点并计算积分
    
17. input=[U(:)';V(:)'];%网络输入
    
18. output=[Tf(:)'];%网络输出
    
19. net=newrb(input,output,0.00001,3*inf_u/nn);%生成逼近网络
    
20. INT=dblquad(@(u,v) sim(net,[u(:)';v*ones(size(u(:)'))]),0,inf_u,0,inf_v);
    
21. 
    
22. function z = Tf_Integrand(u,v)
    
23. m=1;
    
24. n=1;
    
25. a = 1;
    
26. b = 1;
    
27. km = 2*m*pi/a;
    
28. kn = 2*n*pi/b;
    
29. 
    
30. f_rel = @(x,y) (1-cos(km*x)).*(1-cos(kn*y)).*cos(u.*x+v.*y);
    
31. f_img = @(x,y) (1-cos(km*x)).*(1-cos(kn*y)).*sin(u.*x+v.*y);
    
32. tmp_rel = dblquad(f_rel, 0,a, 0,b);
    
33. tmp_img = dblquad(f_img, 0,a, 0,b);
    
34. 
    
35. z = (tmp_rel.^2+tmp_img.^2)./sqrt(u.^2+v.^2);
    
36. end
    
37. end

复制代码

和上面第3楼结果对比：

1. tic,[INT,intt]=alber2008(40,40,30);toc
   
2. Elapsed time is 260.640272 seconds.
   
3. >> INT
   
4. INT =
   
5.     6.9354
   
6. 
   
7. tic,[INTnet,TfInet]=rbf_calc(40,40,30);toc
   
8. NEWRB, neurons = 0, SSE = 0.45778
   
9. NEWRB, neurons = 25, SSE = 0.047086
   
10. NEWRB, neurons = 50, SSE = 0.00730593
    
11. NEWRB, neurons = 75, SSE = 0.00114306
    
12. NEWRB, neurons = 100, SSE = 0.000899809
    
13. NEWRB, neurons = 125, SSE = 0.000185916
    
14. NEWRB, neurons = 150, SSE = 3.66627e-005
    
15. Elapsed time is 121.612190 seconds.
    
16. >> INTnet
    
17. INTnet =
    
18.     7.0517
    
19. 
    
20. tic,INT=rbf_dblquad(40,40,15);INT,toc
    
21. NEWRB, neurons = 0, SSE = 0.138715
    
22. NEWRB, neurons = 25, SSE = 0.00270421
    
23. NEWRB, neurons = 50, SSE = 0.000265526
    
24. NEWRB, neurons = 75, SSE = 4.4852e-005
    
25. INT =
    
26.    6.949356491559914
    
27. Elapsed time is 53.007931 seconds.

复制代码

速度又增快不少。值得说明的是，Tf_Integrand函数在(0,0)点是奇异的，不能选取(0,0)点作为训练样本，限于神经网络外推能力的有限性，(0,0)点附近，网络给出的值存在误差，这些误差主要造成最后的积分误差。

[ 本帖最后由 rocwoods 于 2007-7-30 17:27 编辑 ]

alber2008

呼呼，终于看完所有大虾的方法，我这两天就不一个一个试了（过阵子有空在好好研究研究）
谈谈我的想法
1.
我现在就用taohe 原来写的那段代码进行计算了，换了台性能好点机器，并把所有用到dblquad函数的精度设在了0.001（默认是1.e-6），这样稍微快点了，原来的那个两次嵌套的已经算完了。目前正在用这个方法算一个 三次嵌套的 二重积分（这个更加复杂，第一次的积分函数还是需要分段的：（，所以确实更慢了，没办法，只能一点点算了）

2.
rocwoods大虾的插值法和我原来的想的思路是一样的，就是用最原始的面积×高再叠加的方法计算积分。不过rocwoods通过粗细网格的方式讲这个方法效率改进了一些，但我试了几个系数，效率和1相比没有太大的提高，就没用了，不过还是谢谢，这个思路很好，而且里面的代码比我原来自己写的也好很多，尤其是通过meshgrid和sum函数比
n*n次的循环快很多了。

3.
rocwoods 给出的RBF网络方法是最吸引我的方法，可惜我不懂RBF网络，还需要进一步研究研究。因为理论分析快做完了，开始做实验了，所以最近没空了，过阵子我一定好好研究研究这个方法，可能对我后续的系统辨识有更大的帮助，希望到时候还能与你交流交流哈：）

4.
bainhome 给出的高斯2点积分法，我没有试，但我试过simpsons方法，在第一次求积分时很快，但嵌入在dblquad函数时使用是还是很慢，不知道怎么回事，而且精读不太高，当然simpsons方法本身没有gauss方法高，但我想原理是一样的

5.
我曾经想用Monte Carlo方法来进行计算，这个肯定是最快的了，可惜我从Matlab Center上下的代码似乎用的不对，
我把代码附在后面了，大虾们可以研究研究。

6.
我帮我最近找到的有关多重积分的matlab资料都传上来，给大伙分享一下，估计对大虾是没啥用了，希望对后来者有些帮助。

最后再次感谢各位的帮助！！！！

alber2008

Gauss-Legendre Integration 高斯积分法——对多项式积分很精确

alber2008

Writing Fast MATLAB Code -- 一本教你如何改善matlab执行效率的书，个人认识非常不错：）

alber2008

Using Matlab for Numerical Evaluation of Multiple Integrals

matlab求多重积分的入门教程

alber2008

simpson方法的.m文件

alber2008

蒙特卡洛方法
和
.m文件

topdisk

真是郁闷我看不明白

yiyeguzhou206

请教大虾，不知道这些方法是否适用于比较大的非线性方程组中的复杂的符号函数的二重积分呢？

yiyeguzhou206

我上面说的这个问题是，非线性方程组中的一个对复杂的符号函数进行二重积分的方程（用int(int)方法无法计算下去，由于积分函数含有比较多的自变量）？困扰很久了，希望大虾指点迷津啊！

gentlefu

也不知道rocwoods还来不来这里了 问 一个问题，第四种方法中，不同区域nn取多少合适阿

rocwoods

这个和被积函数的特性有关系，一般来说，变化比较剧烈的nn取的就较大些。nn取的越大，计算时间就越长。实际中，可以让nn从较小值开始，逐渐增大，看相应积分结果的变化。两次积分结果变化较小时或者权衡下所花的时间，来确定nn的值。

shunfly

积分采样点的选取，涉及计算时间和空间，其实，就和自适应高斯积分等等的思路一样
rocwoods兄对震荡函数积分有琢磨么？

这个和被积函数的特性有关系，一般来说，变化比较剧烈的nn取的就较大些。nn取的越大，计算时间就越长。实际中，可以让nn从较小值开始，逐渐增大，看相应积分结果的变化。两次积分结果变化较小时或者权衡下所花的时间 ...
rocwoods 发表于 2009-10-11 14:05

shunfly