#!/usr/bin/env python # ---------------------------------------------------------------- # John Kerl # kerl.john.r@gmail.com # 2006-05-19 # # Contour integral around a (closed) curve C of f(z)/(z-z0)^n dz. # # Example: # # for f in 1 sin cos exp log; do # for p in 0 1 2 3 4; do # echo -n "\$f \$p: " # numintc1 -f \$f -p \$p # done # echo "" # done # ---------------------------------------------------------------- # This software is released under the terms of the GNU GPL. # Please see LICENSE.txt in the same directory as this file. # ================================================================ from __future__ import division # 1/2 = 0.5, not 0. # Import the cmath module rather than the math module so that exp et al. may # take complex arguments. from cmath import * import sys # ---------------------------------------------------------------- # Path def pathz(t, z0): R = 0.01 return z0 + R * exp(1j * t) # ---------------------------------------------------------------- # Non-library-function integrands. def one_func(z): return 1.0 # ---------------------------------------------------------------- def usage(): print >> sys.stderr, "Usage: %s [options]" % (sys.argv[0]) print >> sys.stderr, "Integrates w(z) = f(z)/(z-z0)^n around a small circle centered at z0." print >> sys.stderr, "Options:" print >> sys.stderr, " -t: Tabulate t, z, dz, w, and accumulated sum over the mash." print >> sys.stderr, " -nt {n}: Specify the number of mesh points." print >> sys.stderr, " -z0 {z0}: Specify z0 in 1/(z-z0)^n in the integrand." print >> sys.stderr, " -p {n}: Specify exponent n in 1/(z-z0)^n in the integrand." print >> sys.stderr, " -f {name}: Specify f(z) of integrand: 1, sin, cos, exp, or log." sys.exit(1) # ---------------------------------------------------------------- # Main script # Integrand parameters z0 = 0.0 power = 1 f = one_func # Mesh over path tlo = 0.0 thi = 2.0 * pi nt = 10000 tabulate = 0 # Command-line parsing argc = len(sys.argv) argi = 1 while (argi < argc): if (sys.argv[argi] == '-t'): tabulate = 1 elif (sys.argv[argi] == '-nt'): if ((argc-argi) < 2): usage() argi += 1 nt = int(sys.argv[argi]) elif (sys.argv[argi] == '-p'): if ((argc-argi) < 2): usage() argi += 1 power = int(sys.argv[argi]) elif (sys.argv[argi] == '-z0'): if ((argc-argi) < 2): usage() argi += 1 z0 = complex(sys.argv[argi]) elif (sys.argv[argi] == '-f'): if ((argc-argi) < 2): usage() argi += 1 if (sys.argv[argi] == '1'): f = one_func elif (sys.argv[argi] == 'one'): f = one_func elif (sys.argv[argi] == 'sin'): f = sin elif (sys.argv[argi] == 'cos'): f = cos elif (sys.argv[argi] == 'exp'): f = exp elif (sys.argv[argi] == 'log'): f = log else: usage() else: usage() argi += 1 # Start of numerical approximation to the integral sum = 0.0 t = tlo dt = (thi - tlo) / nt for it in range(0, nt): t1 = t t2 = t + dt z = pathz(t1, z0) dz = pathz(t2, z0) - z w = f(z) / ((z - z0) ** power) sum += w * dz if tabulate: print ">> t = %8.4f/2pi z = %8.4f +%8.4fi dz = %8.4f +%8.4fi w = %8.4f +%8.4fi sum = %8.4f +%8.4fi" % \ (t/2/pi, z.real, z.imag, dz.real, dz.imag, w.real, w.imag, sum.real, sum.imag) t += dt print "%11.7f +%11.7fi" % (sum.real, sum.imag)