myintegrate: Complex integration

In RobinHankin/elliptic: Weierstrass and Jacobi Elliptic Functions

Complex integration

Description

Integration of complex valued functions along the real axis (myintegrate()), along arbitrary paths (integrate.contour()), and following arbitrary straight line segments (integrate.segments()). Also, evaluation of a function at a point using the residue theorem (residue()). A vignette (“residuetheorem”) is provided in the package.

Usage

myintegrate(f, lower,upper, ...)
integrate.contour(f,u,udash, ...)
integrate.segments(f,points, close=TRUE, ...)
residue(f, z0, r, O=z0, ...)

Arguments

f	function, possibly complex valued
lower,upper	Lower and upper limits of integration in myintegrate(); real numbers (for complex values, use integrate.contour() or integrate.segments())
u	Function mapping [0,1] to the contour. For a closed contour, require that u(0)=u(1)
udash	Derivative of u
points	In function integrate.segments(), a vector of complex numbers. Integration will be taken over straight segments joining consecutive elements of points
close	In function integrate.segments(), a Boolean variable with default TRUE meaning to integrate along the segment from points[n] to points[1] in addition to the internal segments
r,O,z0	In function residue() returns f(z0) by integrating f(z)/(z-z0) around a circle of radius r and center O
...	Extra arguments passed to integrate()

Author(s)

Robin K. S. Hankin

Examples

f1 <- function(z){sin(exp(z))}
f2 <- function(z,p){p/z}

myintegrate(f1,2,3)  # that is, along the real axis


integrate.segments(f1,c(1,1i,-1,-1i),close=TRUE)   # should be zero

# (following should be pi*2i; note secondary argument):
integrate.segments(f2,points=c(1,1i,-1,-1i),close=TRUE,p=1)



# To integrate round the unit circle, we need the contour and its
# derivative:

 u <- function(x){exp(pi*2i*x)}
 udash <- function(x){pi*2i*exp(pi*2i*x)}

# Some elementary functions, for practice:

# (following should be 2i*pi; note secondary argument 'p'):
integrate.contour(function(z,p){p/z},u,udash,p=1)      
integrate.contour(function(z){log(z)},u,udash)         # should be -2i*pi
integrate.contour(function(z){sin(z)+1/z^2},u,udash)   # should be zero



# residue() is a convenience wrapper integrating f(z)/(z-z0) along a
# circular contour:

residue(function(z){1/z},2,r=0.1)  # should be 1/2=0.5



# Now, some elliptic functions:
g <- c(3,2+4i)
Zeta <- function(z){zeta(z,g)}
Sigma <- function(z){sigma(z,g)}
WeierstrassP <- function(z){P(z,g)}

jj <- integrate.contour(Zeta,u,udash) 
abs(jj-2*pi*1i)                              # zero to numerical precision
abs(integrate.contour(Sigma,u,udash))        # zero to numerical precision
abs(integrate.contour(WeierstrassP,u,udash)) # zero to numerical precision


# Now integrate f(x) = exp(1i*x)/(1+x^2) from -Inf to +Inf along the
# real axis, using the Residue Theorem.  This tells us that integral of
# f(z) along any closed path is equal to pi*2i times the sum of the
# residues inside it.  Take a semicircular path P from -R to +R along
# the real axis, then following a semicircle in the upper half plane, of
# radius R to close the loop.  Now consider large R.  Then P encloses a
# pole at +1i [there is one at -1i also, but this is outside P, so
# irrelevant here] at which the residue is -1i/2e.  Thus the integral of
# f(z) = 2i*pi*(-1i/2e) = pi/e along P; the contribution from the
# semicircle tends to zero as R tends to infinity; thus the integral
# along the real axis is the whole path integral, or pi/e.

# We can now reproduce this result analytically.  First, choose an R:
R <- 400

# now define P.  First, the semicircle, u1:
u1     <- function(x){R*exp(pi*1i*x)}
u1dash <- function(x){R*pi*1i*exp(pi*1i*x)}

# and now the straight part along the real axis, u2:
u2     <- function(x){R*(2*x-1)}
u2dash <- function(x){R*2}

# Better define the function:
f <- function(z){exp(1i*z)/(1+z^2)}

# OK, now carry out the path integral.  I'll do it explicitly, but note
# that the contribution from the first integral should be small:

answer.approximate <-
    integrate.contour(f,u1,u1dash) +
    integrate.contour(f,u2,u2dash) 

# And compare with the analytical value:
answer.exact <- pi/exp(1)
abs(answer.approximate - answer.exact)


# Now try the same thing but integrating over a triangle, using
# integrate.segments().  Use a path P' with base from -R to +R along the
# real axis, closed by two straight segments, one from +R to 1i*R, the
# other from 1i*R to -R:

abs(integrate.segments(f,c(-R,R,1i*R))- answer.exact)


# Observe how much better one can do by integrating over a big square
# instead:

abs(integrate.segments(f,c(-R,R,R+1i*R, -R+1i*R))- answer.exact)


# Now in the interests of search engine findability, here is an
# application of Cauchy's integral formula, or Cauchy's formula.  I will
# use it to find sin(0.8):

u     <- function(x){exp(pi*2i*x)}
udash <- function(x){pi*2i*exp(pi*2i*x)}

g <- function(z){sin(z)/(z-0.8)}

a <- 1/(2i*pi)*integrate.contour(g,u,udash)


abs(a-sin(0.8))

RobinHankin/elliptic documentation built on Feb. 21, 2024, 6:28 a.m.