轮廓积分C++算法

Algorithm for contour integration C++

本文关键字：算法 C++ 轮廓      更新时间：2023-10-16

我正在尝试编写一个应用数学程序，该程序将计算复平面中的轮廓积分。首先，我想为梯形方法编写一个算法，但我有点难以理解它会是什么样子。毕竟 - 我们通常认为梯形方法与 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;
}

在这里，"函数"只计算 f（z） = 1/z（我确定这是正确完成的），而复笛卡尔是我以 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） 的数量。如果您有任何建议，我将不胜感激。

谢谢。

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

你的梯形规则并没有真正计算复杂的梯形规则，而是真实和复杂之间的一些弗兰肯斯坦。

下面是一个利用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;
}

附言。如果要关心性能，那么integrate会将所有输入作为模板参数。

我喜欢

这里发布的两种解决方案，但我想出的另一个解决方案是参数化我的复杂坐标并在极地工作。由于（在这种情况下）我只在两极周围的小圆圈上求值，因此我的坐标的极坐标形式只有一个变量（θ）。这会将 4D 梯形规则转换为花园品种的 2D 规则。当然，如果我想沿着非圆形轮廓进行积分，这是行不通的，为此我需要上述解决方案。

相关文章：