Algorithm for contour integration C++

Question

I am trying to write an applied math program that will compute contour integrals in the complex plane. For starters, I wanted to write up an algorithm for the trapezoidal method, but I'm somewhat stuck on understand what that would look like. After all - we normally think of the trapezoidal method as for 2D graphs, and here we have f: C -> C so we're talking 4D.

Eventually I'm hoping to compute residues with this algorithm, but when I try the simplest of simple f(z) = 1/z with a contour as a circle of radius 1 around the origin I get nothing near 1 (which is what I should get). Here's my code for the trapezoidal method:

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;
}

Here, 'function' just computes f(z) = 1/z (I'm sure that this is done correctly) and complexCartesian is my class for complex points in the a + bi format:

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();
};

I'm feeling pretty confident that each of these functions is doing what it should be. (I know that there is a standard class for complex numbers but I figured I'd write my own for practice). To actually integrate, I use the following:

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;

When I integrate along the real axis only, or when I replace my function with a constant, I get correct results. Otherwise, I tend to get numbers on the order of 10^(-10) to 10^(-15). If you have any suggestions, I would much appreciate them.

Thanks.

Answer 1

Check your trapezoidal rule: Indeed, contour integrals are defined as a limit of a Riemann sum, $\sum_k f(z_k) \delta z_k$, where the multiplication is to be understood as a complex multiplication, which is not what you do.

Answer 2

Your trapezoidal rule doesn't really compute complex trapezoidal rule, but some Frankenstein between real and complex.

Below is a self-contained example leveraging std::complex, and working "correctly".

#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. If one were to care about performance, then integrate would be taking all the inputs as template parameters.

Answer 3

I like both solutions posted here, but another one that I came up with was to parametrize my complex coordinates and work in polar. Since (in this case) I'm only evaluating at small circles around the poles, the polar form of my coordinates has only one variable (theta). This turns the 4D trapezoidal rule into the garden-variety 2D one. Of course, this doesn't work if I want to integrate along a non-circular contour, for which I would need the above solutions.