轮廓整合算法C++

Question

我想写一个应用数学程序来计算复平面中的轮廓积分。对于初学者，我想写一个梯形方法的算法，但我有点困惑，理解那将是什么样子。毕竟 - 我们通常认为梯形方法与二维图形相同，在这里我们有f：C - > C，所以我们说的是4D。

最终我希望用这个算法计算残基，但是当我尝试简单的f（z）= 1/z的轮廓作为围绕原点的半径为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（我敢肯定，这是正确完成），并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规则。当然，如果我想沿非圆形轮廓进行集成，那么这将不起作用，对此我需要上述解决方案。