计算科学 - 具有无限限制的二维数值积分 (C++) - 吾爱随笔录

具有无限限制的二维数值积分 (C++)

计算科学 C++ 正交一体化

2021-12-11 03:55:11

为了整合形式的二维函数

\int_{1}^{\infty} \int_{- \sqrt{x^{2} - 1}}^{\sqrt{x^{2} - 1}} e^{- x} d y d x,

$\int_{1}^\infty \int_{-\sqrt{x^2-1}}^{\sqrt{x^2-1}} e^{-x} \rm{d}y \rm{d}x,$

我一直在尝试使用以下代码（用 C++ 编写），主要取自 Numerical Recipes 书中，它调用了一个高斯正交例程进行集成：

 static float xsav;
 static float (*nrfunc)(float,float);

 float quad2d(float (*func)(float, float), float x1, float x2)
 {
     float qgaus(float (*func)(float), float a, float b); 
     float f1(float x);

     nrfunc=func;
     return qgaus(f1,x1,x2); 
 }

 float f1(float x) 
 {
     float qgaus(float (*func)(float), float a, float b); 
     float f2(float y);
     float yy1(float),yy2(float);

     xsav=x;
     return qgaus(f2,yy1(x),yy2(x)); 
 }

 float f2(float y)  
 {
     return (*nrfunc)(xsav,y);
 }

此代码适用于具有有限限制的二维积分，但由于外部限制为无穷大而失败。为了解决这个问题，我尝试使用变量的变化：

#define FUNC(x) ((*funk)(-log(x))/(x))

float qgaus(float (*funk)(float), float aa, float bb)
{
    int j;
    float xr,xm,dx,s,a,b;

    b=exp(-aa); 
    a=0.0;

    static float x[]={0.0,0.1488743389,0.4333953941,
        0.6794095682,0.8650633666,0.9739065285};
    static float w[]={0.0,0.2955242247,0.2692667193, 
        0.2190863625,0.1494513491,0.0666713443};

    xm=0.5*(b+a); 
    xr=0.5*(b-a);
    s=0;
    for (j=1;j<=5;j++)
    {
        dx=xr*x[j];
        s += w[j]*(FUNC(xm+dx)+FUNC(xm-dx)); 
    }
    return s *= xr; 
}


float f(float x, float y)
{
    float a = exp(-x);
    return a;
}

float yy1(float x)
{
    float y = -sqrt(x*x-1);
    return y;   
}

float yy2(float x)
{
    float y = sqrt(x*x-1);
    return y;
}



static float xsav;
static float (*nrfunc)(float, float);

float quad2d(float (*func)(float, float), float x1, float x2)
{
    float qgaus(float (*func)(float), float aa, float bb);
    float f1(float x);

    nrfunc=func;
    float t = qgaus(f1,x1,x2);
    return t;
}

float f1(float x)
{
    float qgaus(float (*func)(float), float aa, float bb);
    float f2(float y);
    float yy1(float);
    float yy2(float);

    xsav=x;
    float r = qgaus(f2,yy1(x),yy2(x));
    return r;
}

float f2(float y) 
{
    float k = (*nrfunc)(xsav,y);
    return k;
}

int main ()
{
    float z;
    z = quad2d(f, 1.0, 20.0);
    cout << z << endl;
}

但这仍然没有给出正确的答案。它应该是 $2 \times \rm{BesselK}[1,1] \approx 1.20381$ ，而是给出 2.15501。

任何关于如何修改此代码以解决无限限制的建议将不胜感激！

1个回答

有一种标准方法可以处理积分中的无限限制。这在包的README 中有cubature很好的解释（你应该使用它来进行低维但多维的集成）。

无限间隔

通过变量的变化可以进行无限或半无限区间的积分。最好在一维中说明这一点。

要计算半无限区间上的积分，您可以执行变量的更改 $x=a+t/(1-t)$ ：

$\int_{a}^{\infty} f (x) d x = \int_{0}^{1} f (a + \frac{t}{1 - t}) \frac{1}{(1 - t)^{2}} d t .$ $\int_a^\infty f(x) dx = \int_0^1 f\biggl(a+\frac{t}{1-t}\biggr) \frac{1}{(1-t)^2} dt.$

对于无限区间，您可以执行变量的更改 $x=t/(1 - t^2)$ ：

$\int_{- \infty}^{\infty} f (x) d x = \int_{- 1}^{1} f (\frac{t}{1 - t^{2}}) \frac{1 + t^{2}}{(1 - t^{2})^{2}} d t .$ $\int_{-\infty}^\infty f(x) dx = \int_{-1}^1 f\biggl(\frac{t}{1-t^2}\biggr) \frac{1+t^2}{(1-t^2)^2} dt.$

注意雅可比因子相乘 $f(\dots)$ 在这两个积分中，以及 $t$ 两种情况下的积分是不同的。

在多维中，只需根据需要在每个维度上单独执行变量的这种更改，将被积函数乘以被转换的每个维度的相应雅可比因子。

在您的特定情况下，您可以轻松计算出 $y$ 手动积分，我会使用 GSL 包来计算与现有正交包的积分，它可以处理无限间隔。

使用

\int_{- \sqrt{x^{2} - 1}}^{\sqrt{x^{2} - 1}} \exp (- x) d y = 2 \exp (- x) \sqrt{x^{2} - 1}

$\int_{-\sqrt{x^2-1}}^{\sqrt{x^2-1}} \exp(-x) dy = 2 \exp(-x) \sqrt{x^2-1}$ 然后有一个

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f(double x, void * params)
{
  return 2*exp(-x)*sqrt(x*x-1);
}

int main()
{
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(100);
  double result, error;

  gsl_function F;
  F.function = &f;

  gsl_integration_qagiu(&F, 1, 1e-7, 1e-7, 100, w, &result, &error);

  printf ("result          = %.18f\n", result);
  printf ("estimated error = %.18f\n", error);
  printf ("intervals       = %zu\n", w->size);

  gsl_integration_workspace_free (w);
}

其它你可能感兴趣的问题

上一篇量化传热问题中的非线性程度下一篇多面体内最大的超长方体