3

複素平面の周回積分を計算する応用数学プログラムを書こうとしています。手始めに、台形法のアルゴリズムを書きたかったのですが、それがどのようになるかを理解するのに少し固執しています。結局のところ、私たちは通常、2Dグラフの場合は台形法を考えますが、ここではf：C-> Cであるため、4Dについて話します。

最終的には、このアルゴリズムで残差を計算したいと思っていますが、原点を中心に半径1の円として等高線を使用して、最も単純な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）のオーダーの数字を得る傾向があります。何か提案があれば、よろしくお願いします。

ありがとう。

4

3 に答える 3

3

台形公式は実際には複雑な台形公式を計算しませんが、実数と複素数の間のいくつかのフランケンシュタインを計算します。

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すべての入力をテンプレートパラメータとして使用します。

于 2012-06-21T22:23:20.347 に答える
3

台形公式を確認してください。実際、周回積分はリーマン和の極限として定義されています。$ \ sum_k f（z_k）\ delta z_k $ここで、乗算は複雑な乗算として理解されますが、これはあなたが行うことではありません。

于 2012-06-21T21:43:58.847 に答える
1

ここに投稿された両方のソリューションが気に入っていますが、思いついた別のソリューションは、複雑な座標をパラメーター化し、極座標で作業することでした。(この場合) 極の周りの小さな円でのみ評価しているので、座標の極形式には 1 つの変数 (シータ) しかありません。これにより、4D 台形規則がさまざまな 2D 規則に変わります。もちろん、上記のソリューションが必要な非円形の輪郭に沿って統合したい場合、これは機能しません。

于 2012-06-22T17:19:21.143 に答える