輪廓整合算法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規則。當然，如果我想沿非圓形輪廓進行集成，那麼這將不起作用，對此我需要上述解決方案。