轮廓积分算法 C++

Question

我正在尝试编写一个应用数学程序来计算复平面中的轮廓积分。对于初学者，我想为梯形方法编写一个算法，但我有点难以理解它会是什么样子。毕竟 - 我们通常将梯形方法视为 2D 图，这里我们有 f: C -> C 所以我们说的是 4D。

最终我希望用这个算法计算残差，但是当我尝试最简单的简单 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) 的数字。如果您有任何建议，我将不胜感激。

谢谢。

Answer 1

检查你的梯形规则：确实，轮廓积分被定义为黎曼和的极限，$\sum_k f(z_k) \delta z_k$，其中乘法被理解为复数乘法，这不是你做。

Answer 2

您的梯形规则并不是真正计算复杂梯形规则，而是计算真实和复杂之间的一些科学怪人。

下面是一个利用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 会将所有输入作为模板参数。

Answer 3

我喜欢这里发布的两种解决方案，但我想出的另一种解决方案是将我的复杂坐标参数化并在极地工作。因为（在这种情况下）我只在极点周围的小圆圈上进行评估，所以我的坐标的极坐标形式只有一个变量（theta）。这将 4D 梯形规则变成了普通的 2D 规则。当然，如果我想沿非圆形轮廓进行积分，这不起作用，我需要上述解决方案。