平面离散点上求二重积分求助

xmzhy500234

我尽量简洁说明
已知了一系列离散点（平均分布）的函数值FX(I,J)，即X,Y平面上某一函数值，区域X：1-M，Y:1-N.
现在要求二重积分
{F}_{x}=\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}_{S}{f}_{x}dxdy
求大神指导下，我粘贴下我自己的程序，其中H是节点距离。

[Fortran] 纯文本查看 复制代码

FEX=0.
   do i=1,M         !!外循环对x方向积分
           FSX(i)=0.
           FSY(i)=0.
           do j=1,N    !!内循环对y方向积分 
              FSX(I)=FSX(i)+(FX(i,j)+FX(i,j+1))*H/2.0
       end do 
    end do
    do i=1,M 
           FEX=FEX+(FSX(i)+FSX(i+1))*H/2.0
    end do

糖盒love玲珑

采样点很稀疏的话, 这个积分值就很难确定, 100 个人可以算出 101 个结果(我可以算2种结果). 在你给的采样点外其他的函数值都是未知的, 该积分无解.

kerb

这个你最好手工先推到以下
辛普森公式是对于区间[a,b],积分
S=\frac{b-a}{6}(f(a)+f(\frac{a+b}{2})+f(b))
如果是等分区间
[a=x_1,x_2,\cdots,x_n=b]，
对区间
[x_1,x_2,x_3], [x_3,x_4,x_5], \cdots, [x_{n-4},x_{n-3},x_{n-2}], [x_{n-2},x_{n-1},x_n]
分别使用辛普森公式，然后加起来，可以得到
\frac{h}{6}(f(x_1)+4f(x_2)+2f(x_3)+4f(x_4)+\cdots+4f(x_{n-1})+f(x_n))
对于二重积分
S=\int_a^b\int_c^d f_2dydx=\int_a^b f_1 dx; f_1=\int_c^d f_2 dy，
其中
\int_a^b f_1=\frac{h_x}{6}(f_1(x_1)+4f_1(x_2)+2f_1(x_3)+4f_1(x_4)+\cdots+4f_1(x_{n-1})+f_1(x_n)),
你把f_1对应的公式写进去，
\int_c^d f_2=\frac{h_y}{6}(f_2(y_1)+4f_2(y_2)+2f_2(y_3)+4f_2(y_4)+\cdots+4f_2(y_{n-1})+f_2(y_n))
然后整理一下，再编程序，另外，上面的复合辛普森公式，你再查一下书，我这里只是顺手写的，似乎还有更好的形式

pasuka

xmzhy500234 发表于 2015-11-23 15:30
问题到重点了
被积函数，主要就是没有函数！，是一系列值，最关键就在这里了。其实是函数就简单多了 ...

RBF径向基拟合，然后数值积分
或者平面划分三角形网格，有背景网格后，再数值积分
可以往无网格方向的文献走走，应该有不少办法

xmzhy500234

pasuka 发表于 2015-11-23 15:14
误差大不大，被积函数的性质很重要，振荡函数、间断函数都有对应的特殊处理
离散点是平均分布，那么蒙特 ...

问题到重点了
被积函数，主要就是没有函数！，是一系列值，最关键就在这里了。其实是函数就简单多了，两分钟就做好。
我这样写可能更明白点，数据格式如下：
X    Y      fx      fy（X坐标，Y坐标，X方向fx，Y方向fy）
...
...
数据下面是循环，比如一个X，对应63个Y，然后X逐步按节点前进。意思就是整个平面的网格点。

pasuka

xmzhy500234 发表于 2015-11-23 14:28
我目前用同一个程序，计算了两个例子，一个例子是63X63节点，一个是127X127，两者应该在结果上是相差不大 ...

误差大不大，被积函数的性质很重要，振荡函数、间断函数都有对应的特殊处理
离散点是平均分布，那么蒙特卡洛方法最合适不过，网格逐步加密肯定能逼近真实值
办法总比问题多，关键得自己动手实践

xmzhy500234

kerb 发表于 2015-11-23 01:09
这个有多种方法，根据你对积分精度要求和区域大小而定，从你的题目上看，你是等间距的，二次一维线积分即可 ...

补充说明：127X127的源数据是在63X63基础上插值得来的。

xmzhy500234

kerb 发表于 2015-11-23 01:09
这个有多种方法，根据你对积分精度要求和区域大小而定，从你的题目上看，你是等间距的，二次一维线积分即可 ...

我目前用同一个程序，计算了两个例子，一个例子是63X63节点，一个是127X127，两者应该在结果上是相差不大的，但我计算的结果是相差很大的

xmzhy500234

pasuka 发表于 2015-11-22 18:49
最简单直白的就是蒙特卡罗方法，必要时改成并行或者CUDA加速也是挺容易的
希望档次高些，可以用Clenshaw-Cu ...

不是我不想提高精度，而是我先用最简单的尝试，算出来结果不对，我始终认为我程序里这个积分公式是不正确的，然而也想不出来该如何做。

kerb

这个有多种方法，根据你对积分精度要求和区域大小而定，从你的题目上看，你是等间距的，二次一维线积分即可，如果要求较高精度线积分可选用矩形、梯形公式、辛普森公式，高阶牛顿-柯兹公式，你的代码用的是体形公式