C++轮廓积分算法

Question

我正在尝试编写一个应用数学程序，用于计算复平面上的等高线积分。首先，我想编写梯形法的算法，但我对其外观有些困惑。毕竟 - 我们通常认为梯形法是针对二维图像的，而这里我们有f：C->C，所以我们在谈论四维。

最终，我希望使用此算法计算残留值，但当我尝试最简单的f（z）= 1 / z并将轮廓线设置为围绕原点的半径为1的圆时，我得不到接近1的值（这是我应该得到的）。以下是我的梯形法代码：

complexCartesian trapezoid(complexCartesian *c1, complexCartesian *c2)
{
     complexCartesian difference = *c1 - *c2;
     complexCartesian ans(0.5 * (function(c1).real + function(c2).real) * 
                                                     difference.mag(), 
                          0.5 * (function(c1).imag + function(c2).imag) *
                                                     difference.mag());
     return ans;
}

这里的 'function' 函数计算的是 f(z) = 1/z (我确定计算正确)，而 complexCartesian 则是用于表示 a + bi 格式下的复数点的类：

class complexCartesian
{
    public:
    double real;
    double imag;

    complexCartesian operator+ (const complexCartesian& c) const;
    complexCartesian operator- (const complexCartesian& c) const;
    complexCartesian(double realLocal, double imagLocal);
    double mag(); //magnitude
    string toString();
    complexPolar toPolar();
};

我很有信心这些函数都能够正常运行。(我知道有一个标准的复数类，但为了练习，我决定自己写。) 要进行积分，我使用以下方法：

double increment = .00001;
double radius = 1.0;
complexCartesian res(0,0); //result
complexCartesian previous(1, 0); //start the point 1 + 0i

for (double i = increment; i <= 2*PI; i+=increment)
{
    counter++;
    complex cur(radius * cos(i), radius * sin(i));
    res = res + trapezoid(&cur, &previous);
    previous = cur;
}

cout << "computed at " << counter << " values " << endl;
cout << "the integral evaluates to " << res.toString() << endl;

当我只沿着实轴进行积分，或者将我的函数替换为常数时，我得到正确的结果。否则，我往往会得到10^(-10)到10^(-15)数量级的数字。如果您有任何建议，我会非常感激。

谢谢。

- alexvas

---

我已经有一段时间没有做过等高线积分了，但是以上的积分是否可能计算为0呢？ - templatetypedef

No. 1/z在原点处有一个极点，且该极点的残留为1（lim_{z\to 0} z(1/z) = 1）。因此，根据残留定理，积分应该计算得出2(pi)(i)。 - alexvas

啊，好的。谢谢你让我知道！我很好奇这是否只是数值不稳定性问题。 - templatetypedef

你不想使用极坐标的特定原因是什么？ - ev-br

Answer 1

你的梯形法并不能真正计算“复杂”的梯形法，而是介于实数和复数之间的一种混合体。

以下是一个使用std::complex的自包含示例，并且可以正确工作。

#include <cmath>
#include <cassert>
#include <complex>
#include <iostream>

using std::cout;
using std::endl;
typedef std::complex<double> cplx;

typedef cplx (*function)(const cplx &);
typedef cplx (*path)(double);
typedef cplx (*rule)(function, const cplx &, const cplx &);

cplx inv(const cplx &z)
{
    return cplx(1,0)/z;
}

cplx unit_circle(double t)
{
    const double r = 1.0;
    const double c = 2*M_PI;
    return cplx(r*cos(c*t), r*sin(c*t));
}

cplx imag_exp_line_pi(double t)
{
    return exp(cplx(0, M_PI*t));
}

cplx trapezoid(function f, const cplx &a, const cplx &b)
{
    return 0.5 * (b-a) * (f(a)+f(b));
}

cplx integrate(function f, path path_, rule rule_)
{
    int counter = 0;
    double increment = .0001;
    cplx integral(0,0);
    cplx prev_point = path_(0.0);
    for (double i = increment; i <= 1.0; i += increment) {
        const cplx point = path_(i);
        integral += rule_(f, prev_point, point);
        prev_point = point;
        counter ++;
    }

    cout << "computed at " << counter << " values " << endl;
    cout << "the integral evaluates to " << integral << endl;
    return integral;
}

int main(int, char **)
{
    const double eps = 1E-7;
    cplx a, b;
    // Test in Octave using inverse and an exponential of a line
    // z = exp(1i*pi*(0:100)/100);
    // trapz(z, 1./z)
    // ans =
    //   0.0000 + 3.1411i
    a = integrate(inv, imag_exp_line_pi, trapezoid);
    b = cplx(0,M_PI);
    assert(abs(a-b) < eps*abs(b));

    // expected to be 2*PI*i
    a = integrate(inv, unit_circle, trapezoid);
    b = cplx(0,2*M_PI);
    assert(abs(a-b) < eps*abs(b));

    return 0;
}

PS. 如果关心性能，那么integrate将所有输入作为模板参数。

- Kuba hasn't forgotten Monica

Answer 2

检查一下你的梯形公式：事实上，轮廓积分被定义为一个黎曼和的极限，即 $\sum_k f(z_k) \delta z_k$，其中乘法应理解为复数乘法，而这不是你所做的。

- ev-br

---

我认为delta z是一个大小。因此使用difference.mag()。你的意思是我应该乘以复数差异而不是实数差异吗？ - alexvas

是的。例如，请参见此处的示例6.6：http://math.fullerton.edu/mathews/c2003/ContourIntegralMod.html - ev-br

是的 - 使用复数进行“计算”的规则是您只需使用复数计算所有内容。 您不能随意地从复数导出实数并任意地计算公式的某些部分。 这就是使用 mag（） 的原因：您在复数中有一个完美的公式，然后选择毫无理由地破坏它 :) - Kuba hasn't forgotten Monica

Answer 3

我喜欢这里发布的两种解决方案，但我想到了另一种方法，那就是将我的复杂坐标参数化并在极坐标下工作。由于（在这种情况下）我只在极点周围的小圆上进行评估，因此我的坐标的极坐标形式只有一个变量（θ）。这将4D梯形法转换为普通的2D梯形法。当然，如果我想沿着非圆形轮廓积分，则无法使用此方法，而需要使用上述解决方案。

- alexvas