"Real Contour Integration": Code Followup

Question

Adrian Keister had a wonderful answer to my question. In short, his answer allowed me to make a form of contour integration that picks up (some of the time) only real roots of a function $f(z)$. I started to code this in Python 3.x, but while I was debugging it, sometimes it would give me a negative answer or that there is no roots in an interval I know a root lies in. I can't see anything wrong with my code, so I'm hoping that someone else will point out my folly.

import cmath as cm
import math as m
start=150
end=181
subdiv=100
front=1/(2*cm.pi*1j)
d=1000
eps=1/d
def f(x):
  return x**2 - 358*x + 32040
def f_1(x):
   return 2*(-179 + x)
def g_1_1(x):
  return start+(x*(end-start))+(1j*eps)
def g_1_2(x):
  return (end-start)

def g_2_1(x):
  return end+(1j*(eps-(2*x*eps)))
def g_2_2(x):
  return -2*1j*eps

def g_3_1(x):
  return end+(x*(start-end))-(1j*eps)
def g_3_2(x):
  return (start-end)

def g_4_1(x):
  return start+(1j*((-1*eps)+(2*x*eps)))
def g_4_2(x):
  return 2*1j*eps

def c(g,g_1,x):
  return (f_1(g(x))/f(g(x)))*g_1(x)

def quad(g,g_1):
  total=0
  for i in range(0,subdiv):
    a=(i)*((1)/subdiv)
    b=(i+1)*((1)/subdiv)
    total+=(b-a)*((c(g,g_1,a)+c(g,g_1,b))/2)
  return total

def g_c():
  total=0
  total+=quad(g_1_1,g_1_2)
  total+=quad(g_2_1,g_2_2)
  total+=quad(g_3_1,g_3_2)
  total+=quad(g_4_1,g_4_2)
  return front*total

print(g_c())

In an interval ($150$ to $190$) where two roots lie very close to each other ( both close to$179$), Python returns 0.0031720503705679384+5.624689839008629e-29j when I thought I would get 1.999986474+0.6746345e-99j

EDIT: Here is the new code:

def g_1_1(x):
  return start+((1j*eps)*(1-(2*x)))
def g_1_2(x):
  return -2*1j*eps

def g_2_1(x):
  return start+(x*(end-start))-(1j*eps)
def g_2_2(x):
  return (end-start)

def g_3_1(x):
  return end+((1j*eps)*(-1+(2*x)))
def g_3_2(x):
  return 2*1j*eps

def g_4_1(x):
  return end+(x*(start-end))+(1j*eps)
def g_4_2(x):
  return (start-end)

def c(g,g_1,x):
  return (f_1(g(x))/f(g(x)))*g_1(x)

def quad(g,g_1):
  total=0
  for i in range(0,subdiv):
    a=(i)*((1)/subdiv)
    b=(i+1)*((1)/subdiv)
    total+=(b-a)*((c(g,g_1,a)+c(g,g_1,b))/2)
  return total

def g_c():
  total=0
  total+=quad(g_1_1,g_1_2)
  total+=quad(g_2_1,g_2_2)
  total+=quad(g_3_1,g_3_2)
  total+=quad(g_4_1,g_4_2)
  return front*total

print(g_c())

Now it gives me an answer of -0.0031720503705679384-0j.

Answer 1

So let's apply the argument principle to the function $f(x)=x^2-358x+32040.$ This function has two known zeros, at $178$ and $180$. So, referring to my other answer, let us take the interval from $150$ to $190$. We would have $a=150,\; b=190.$ Let's use the circular contour, since we know there are no poles. We would have \begin{align*} f(z)&=z^2-358z+32040\\ f'(z)&=2z-358\\ c&=170\\ r&=20\\ z&=170+20e^{i\theta}\\ dz&=20\,i\,e^{i\theta}\,d\theta. \end{align*} The necessary contour integral becomes \begin{align*} Z-P&=\frac{1}{2\pi i}\oint_C\frac{f'(z)}{f(z)}\,dz\\ &=\frac{1}{2\pi i}\oint_C\frac{2z-358}{z^2-358z+32040}\,dz\\ &=\frac{1}{2\pi i}\int_0^{2\pi}\frac{2(170+20e^{i\theta})-358}{(170+20e^{i\theta})^2-358(170+20e^{i\theta})+32040}\,20\,i\,e^{i\theta}\,d\theta\\ &=\frac{10}{\pi}\int_0^{2\pi}\frac{40e^{i\theta}-18}{80-360e^{i\theta}+400e^{2i\theta}}\,e^{i\theta}\,d\theta\\ &=\frac{10}{\pi}\int_0^{2\pi}\frac{20e^{i\theta}-9}{40e^{-i\theta}-180+200e^{i\theta}}\,d\theta\\ &=\frac{10}{\pi}\cdot\frac{\pi}{5}\\ &=2, \end{align*} as required! Similarly, if we take the rectangle contour \begin{array}{|c|c|c|c|} \hline &z &t\;\text{interval} &dz\\ \hline \gamma_1 &a+i\varepsilon(1-2t) &[0,1] &-2i\varepsilon\,dt \\ \hline \gamma_2 &a+t(b-a)-i\varepsilon &[0,1] &(b-a)\,dt \\ \hline \gamma_3 &b+i\varepsilon(-1+2t) &[0,1] &2i\varepsilon\,dt \\ \hline \gamma_4 &b+t(a-b)+i\varepsilon &[0,1] &(a-b)\,dt \\ \hline \end{array} this translates to \begin{array}{|c|c|c|c|} \hline &z &t\;\text{interval} &dz\\ \hline \gamma_1 &150+i\varepsilon(1-2t) &[0,1] &-2i\varepsilon\,dt \\ \hline \gamma_2 &150+40t-i\varepsilon &[0,1] &40\,dt \\ \hline \gamma_3 &190+i\varepsilon(-1+2t) &[0,1] &2i\varepsilon\,dt \\ \hline \gamma_4 &190-40t+i\varepsilon &[0,1] &-40\,dt \\ \hline \end{array} For an example, I will take the $\gamma_1$ integral. Let $\gamma=\gamma_1\cup\gamma_2\cup\gamma_3\cup\gamma_4.$ For reference: \begin{align*} Z-P &=\frac{1}{2\pi i}\oint_\gamma\frac{f'(z)}{f(z)}\,dz \\ &=\frac{1}{2\pi i}\left[\sum_{k=1}^4\int_{\gamma_k}\frac{f'(z)}{f(z)}\,dz\right]. \end{align*} So we have \begin{align*} \int_{\gamma_1}\frac{f'(z)}{f(z)}\,dz &=-2i\,\varepsilon\int_0^1\frac{f'(150+i\varepsilon(1-2t))}{f(150+i\varepsilon(1-2t))}\,dt\\ &=-2i\,\varepsilon\int_0^1\frac{2(150+i\varepsilon(1-2t))-358}{(150+i\varepsilon(1-2t))^2-358(150+i\varepsilon(1-2t))+32040}\,dt\\ &=\ln\left(\varepsilon^2-58i\varepsilon-840\right)-\ln((\varepsilon+28i)(\varepsilon+30i))+2i\pi. \end{align*} Here are the results for the other segments: \begin{align*} \int_{\gamma_2}\frac{f'(z)}{f(z)}\,dz&=\ln[(\varepsilon+10i)(\varepsilon+12i)]-\ln[\varepsilon^2-58\varepsilon i-840]\\ \int_{\gamma_3}\frac{f'(z)}{f(z)}\,dz&=\ln(\varepsilon^2-22\varepsilon i-120)-\ln[(\varepsilon+10i)(\varepsilon+12i)]+2\pi i\\ \int_{\gamma_4}\frac{f'(z)}{f(z)}\,dz&=\ln[(\varepsilon+28i)(\varepsilon+30i)]-\ln(\varepsilon^2-22\varepsilon i-120). \end{align*} This makes the sum equal to $4\pi i,$ (all the logarithm terms cancel out) which yields the correct result when you divide by $2\pi i.$

You can use these individual results to debug your code.