Contour integration in Racket

I created the first half of this last year for pi day. Back then, I wrote it in C and just integrated `1/z` along a square contour to get π. A few days ago I came back to it, and rewrote it in racket. But then I realized, I can actually make a function that integrates along any curve very simply.

---

## Polygonal contour

I want to write a function, that you can pass a function and a list of point, which will integrate the function along a polygonal path defined by the points passed into the function.

First, I want to define a function, that takes two points and integrates along the line defined by these points.
I will just use the rectangle method.

I will simply divide the contour into `n` pieces and treat the function on each piece as constant. Provided continuity, we can assure that as we take the number of pieces to the limit, we will converge to an answer. After all, this is simply the Riemann definition of the integral.

This leads us to the following function definition:


    (define (integrate f a b steps)
      (define dx (/ (- b a) steps))
      (for/sum ([i steps])
        (* dx (f (+ a (* dx i))))))

In this function definition, I'm defining a function `integrate`. This function takes a function, two point and a number of steps.
From there, we define the length of each step.

    (define dx (/ (- b a) steps)

or in infix notation

    dx = (b-a)/steps

Then, I perfom a for/sum. Simply said, it's a for loop, where the output of each call inside the loop (or more precisely, each _last_ result of the loop body) is accumulated. The result of this accumulation is the return value of the loop.

Inside the loop,

    (* dx (f (+ a (* dx i))))

or in infix notation

    dx * f(a + (i*dx))

we take the "height" of our rectangle and multiply it with its width `dx`. The loop goes over the list 0-steps.

## Setting up the contour

We want to create a function, that takes our function to integrate, our selected number of steps, and our list of points.


    (require racket/list)
    (define (contour f steps . points)
      (define usable-points (append points (list (first points))))
      (for/sum ([i (sub1 (length usable-points))])
            (integrate f (list-ref usable-points i) (list-ref usable-points (add1 i)) steps)))

In this snippet, we are adding a requirement for the `racket/list` library, to allow us to use functions `first` and `last`.

We are doing a very similar thing to the snippet above, but before we can iterate over our list of points, we have to finagle it a bit.
We simply add our first point back to the end of the list.

This way, the last step of the integration closes the loop, and makes our contour into a closed curve.

In detail:

We are defining a function `contour` which takes a function to be integrated, a number of steps for each line of the polygon, and a variable length list of points defining our contour.
Then, we are again doing a `for/sum` and each time calling the function `integrate` we defined earlier and feeding it a pair of point from the list supplied to `contour`. Here, we can see the finagling click into place. On each evaluation of the body of the `for/sum` loop, we are feeding it the point pairwise, and to close the loop, we have to make the last evaluation be with the last and first point in that order.

## Putting it all together

We end up with the following code all together:


    (require racket/list)
    
    (define (integrate f a b steps)
      (define dx (/ (- b a) steps))
      (for/sum ([i steps])
        (* dx (f (+ a (* dx i))))))
    
    (define (contour f steps . points)
      (define usable-points (append points (list (first points))))
      (for/sum ([i (sub1 (length usable-points))])
            (integrate f (list-ref usable-points i) (list-ref usable-points (add1 i)) steps)))

I moved the `required` call to the top of the code, but otherwise, this is a simple copy-paste of the cells in the order as they were presented.

Now to see it in action, we can integrate the function `1/x` on the "unit square" i.e. the square defined by the points `1 + i`, `-1 + i`, `-1 - i` and `-1 + i`. First, we define the function to be integrated:


    (define (oneover x)
      (/ 1 x))

And we call the `contour` function with the points described above:


    (contour oneover 100 1.+i -1.+i -1.-i 1.-i)




`-0.03999999999999962+6.283051973846504i`



I called the function with 100 steps on each line segment. The expected value of this integral is `2πi` and we see we are approaching this value. I'll clean up the output by isolating the imaginary part and dividing by two, giving us an expected value of just π. Let's try more steps.


    (/ (imag-part (contour oneover 1000 1.+i -1.+i -1.-i 1.-i)) 2)



`3.1415919869231255`




    (/ (imag-part (contour oneover 5000 1.+i -1.+i -1.-i 1.-i)) 2)




`3.1415926269231225`




    (/ (imag-part (contour oneover 10000 1.+i -1.+i -1.-i 1.-i)) 2)




`3.1415926469231286`



Which is pretty damn close to π

# Arbitrary contour

First, we can actually remove all the finagling of the point array by simply requiring, that if the contour should be closed, the first point should also be added at the end, this does put extra burden on the usage of the function, but it actually allows us to use the exact same function for open and closed contours


    (require racket/list)
    
    (define (integrate f a b steps)
      (define dx (/ (- b a) steps))
      (for/sum ([i steps])
        (* dx (f (+ a (* dx i))))))
    
    (define (contour f steps . points)
      (for/sum ([i (sub1 (length points))])
            (integrate f (list-ref points i) (list-ref points (add1 i)) steps)))
    ;; here I'm defining the function 1/x and integrating along the square contour as before
    (define (f x) (/ 1 x))
    (/ (imag-part (contour f 10000 1.+i -1.+i -1.-i 1.-i 1.+i)) 2) ; notice the first point is also the last point




`3.1415926469231286`



## Line integrals done properly

In the previous examples we were dividing the line directly and manipulating everything like that.

But we can define a function γ which will be a function from some parameter space to $\mathbb{C}$, which will basically parametrize the given curve. Then, when we're building the Riemann sum, we aren't partitioning the curve directly anymore, we're actually partitioning the parameter space (if thinking about this physically, the function γ could describe a position of a particle at a certain time, so we aren't partitioning the position of the particle anymore, we're partitioning the time parameter) so we need to then carry this partitioning from the parameter space forward to the curve described by γ.

You could think of this like a Riemann-Stieltjes integral I guess.

So the integral has the form

![contour integral as a Riemann sum](https://tilde.team/~fanoplanes/blog/equation.svg)

I am creating a function `arbitrary-integrate` that takes my function to be integrated, my curve parametrization, the start and end point of the parameter, and my number of steps.


    (require racket/math)
    (define (arbitrary-integrate f gamma start stop steps)
      (define dt (/ (- stop start) steps))
      (for/sum ([i steps])
        (* (- (gamma (+ start (* dt (add1 i)))) 
              (gamma (+ start (* dt i)))) 
           (f (gamma (+ start (* dt i)))))))
    
    (define (circle t) ; here I define a circular contour, the parameter ranges from 0 to 2*pi
      (exp (* 0+i t)))
    (define (f z) ; I will be again integrating 1/z
      (/ 1 z))
    (/ (imag-part (arbitrary-integrate f circle 0 (* 2 pi) 10000)) ; and I am extracting the imaginary part and dividing by 2 to get pi
       2)




`3.1415924468808663`



So the whole thing is just


    (require racket/math)
    
    (define (arbitrary-integrate f gamma start stop steps)
      (define dt (/ (- stop start) steps))
      (for/sum ([i steps])
        (* (- (gamma (+ start (* dt (add1 i)))) 
              (gamma (+ start (* dt i)))) 
           (f (gamma (+ start (* dt i)))))))