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

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

xmzhy500234

FEX是最终结果，FSX(I)是中间变量

pasuka

最简单直白的就是蒙特卡罗方法，必要时改成并行或者CUDA加速也是挺容易的
希望档次高些，可以用Clenshaw-Curtis积分，同样可以借助FFTW加速计算
或者Gauss-Kronrod积分，对于振荡函数挺管用的

kerb

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

xmzhy500234

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

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

xmzhy500234

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

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

xmzhy500234

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

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

pasuka

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

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

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 15:30
问题到重点了
被积函数，主要就是没有函数！，是一系列值，最关键就在这里了。其实是函数就简单多了 ...

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