등고선 통합을위한 알고리즘 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; 
} 

는 '기능'은 바로 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

나는 두 가지 솔루션이 여기에 게시되는 것을 좋아하지만, 또 다른 하나는 필자의 복잡한 좌표를 매개 변수화하고 극지방에서 작업하는 것입니다. (이 경우) 극점 주위의 작은 원에서만 평가하기 때문에 극좌표 형태의 좌표는 하나의 변수 (세타)를 갖습니다. 그러면 4D 사다리꼴 규칙이 정원 종류 2D로 바뀝니다. 물론, 위의 솔루션을 필요로하는 비 원형 윤곽을 따라 통합하려는 경우에는 이것이 작동하지 않습니다.