myintegrate.Rd
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.
myintegrate(f, lower,upper, ...)
integrate.contour(f,u,udash, ...)
integrate.segments(f,points, close=TRUE, ...)
residue(f, z0, r, O=z0, ...)
function, possibly complex valued
Lower and upper limits of integration in myintegrate()
;
real numbers (for complex values, use integrate.contour()
or
integrate.segments()
)
Function mapping \([0,1]\) to the contour. For a closed contour, require that \(u(0)=u(1)\)
Derivative of u
In function integrate.segments()
, a vector of complex
numbers. Integration will be taken over straight segments joining
consecutive elements of points
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
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()
f1 <- function(z){sin(exp(z))}
f2 <- function(z,p){p/z}
myintegrate(f1,2,3) # that is, along the real axis
#> [1] 0.0551893+0i
integrate.segments(f1,c(1,1i,-1,-1i),close=TRUE) # should be zero
#> [1] 0+0i
# (following should be pi*2i; note secondary argument):
integrate.segments(f2,points=c(1,1i,-1,-1i),close=TRUE,p=1)
#> [1] -2.596731e-18+6.283185i
# 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)
#> [1] -3.561641e-17+6.283185i
integrate.contour(function(z){log(z)},u,udash) # should be -2i*pi
#> [1] -8.881784e-16-6.283185i
integrate.contour(function(z){sin(z)+1/z^2},u,udash) # should be zero
#> [1] 1.079679e-16-8.298917e-15i
# 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
#> [1] 0.5-7.160673e-17i
# 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
#> [1] 9.302781e-16
abs(integrate.contour(Sigma,u,udash)) # zero to numerical precision
#> [1] 2.576717e-12
abs(integrate.contour(WeierstrassP,u,udash)) # zero to numerical precision
#> [1] 2.482534e-16
# 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)
#> [1] 6.244969e-07
# 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)
#> [1] 5.157772e-07
# 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)
#> [1] 2.319341e-08
# 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))
#> [1] 3.752575e-14