Summation in the form of Jacobi Theta Function

elliptic functionspythonsummationtheta-functions

TL;DR: My summation should give the same result when it is expressed as the Jacobi Theta function. It gives the same results for some set of inputs but then gives exactly the half of it for other set of inputs. Where do I go wrong?

I have a summation of the form

\begin{equation}
\tag{1}
\sum\limits_{i=0}^\infty \frac{1}{A^x}\left[\exp\left(-\frac{(id-a)^2}{A}\right)+\exp\left(-\frac{(id+a)^2}{A}\right)\right]
\end{equation}

So, I proceed as

\begin{align}
=&\frac{1}{A^x}\sum\limits_{i=0}^\infty \exp\left(-\frac{i^2d^2-2iad+a^2}{A}\right)\left(\exp\left(-\frac{i^2d^2+2aid+a^2}{A}\right)\right)\\
=&\frac{1}{A^x}\sum\limits_{i=0}^\infty \exp\left(-\frac{i^2d^2}{A}\right)\exp\left(-\frac{a^2}{A}\right)\left(\exp\left(-\frac{2aid}{A}\right)+\exp\left(\frac{2aid}{A}\right)\right)\\
\tag{2}
=&\frac{2}{A^x}\exp\left(-\frac{a^2}{A}\right)\sum\limits_{i=0}^\infty \exp\left(-\frac{i^2d^2}{A}\right)\cosh\left(\frac{2aid}{A}\right),
\end{align}
where I used $2\cosh(x)=\exp(x)+\exp(-x)$.

As a side note, $\cosh$ substitution doesn't help numerically because it goes too high too fast while $\exp$ term goes too low too fast so some precision might be lost. But with the $\cosh$ substitution, this looks more like the third Jacobi Theta function, i.e.,

$$
\tag{3}
=\frac{2}{A^x}\exp\left(-\frac{a^2}{A}\right)\vartheta_3\left(\exp\left(-\frac{d^2}{A}\right),\frac{ad}{A}j\right),
$$
where $j=\sqrt{-1}$ and $\vartheta_3(q,z)=1+2\sum_{i=1}^\infty q^{i^2}\cos(2iz)$, as it is defined here. Using (3) instead of (1) is not great as well, because in spite of how fancy it looks, as $q\rightarrow 1$, $\vartheta_3$ calculations go haywire, while choosing an arbitrarily large $N$ to replace $\infty$ in (1) solves the issue even when $q=1-\epsilon$.

And the question is here: why is it that (1) and (3) do not give me the same value for every $A$?

Things that I eliminated so far:

  • Too small $N$ for $\infty$: I used both $N=10000$ and $N=50000$ in the plots, and since they are too aligned, $N=10000$ line is not visible.
  • mpmath.dps too low: I set mpmath.dps=50000, and it is still unchanged.
  • mpmath.jtheta might be inaccurate as $q\rightarrow 1$: Wolfram alpha gives the same results as mpmath verbatim.
  • Maybe $a$ plays some tricks: If $a=0$, we would have $z=0$ and $\vartheta_3(\exp(-d^2/A),0)$, i.e., $\cos$ term in the Jacobi Theta function would be eliminated. But even with $a=0$, the results are unchanged.
  • Maybe $x$ is the problem: Not really, I tried $x \in [0,0.1,0.5,1,2]$ and still no change. The results are identical.

My hunches that I cannot test myself:

  • I make an implicit assumption that becomes visible only when $A$ is large enough.
  • It is about machine precision that makes both wolfram alpha and mpmath act weird after a certain point, especially if they use the same algorithm to compute $\vartheta_3(.)$. I tried to check which algorithm they use but it was not readily available.

Also, the result with (1) is exactly the half of (3), which is also weird and works against my second hunch.

So, where do I go wrong?

Appendix

Here, the figures are $A$ vs the result of (1) divided by (3) for the same input set. $A$ values are given in the x-axis.

Plots of (1)/(3) for $a=0.1$ for different $x$ values

Plots of (1)/(3) for <span class=$a=0.1$ for different $x$" />

Plots of (1)/(3) for $a=0$ for different $x$ values

Plots of (1)/(3) for <span class=$a=0$ for different $x$" />

d=1
N=10000 #also 50000
a=0.1 #also 0
AA=np.logspace(-4,15,100) #the array that carries A values
pow=0 #[0,0.1,0.5,1,2]
def summation(A,d,a,N,pow):
    s=0
    for i in range(N):
        b1=i*d-a
        b2=i*d+a
        s+=np.exp(-(b1**2)/A)/(A**pow)
        s+=np.exp(-(b2**2)/A)/(A**pow)
    return s

def vartheta(A,d,a,pow):
    p=(a*d)/A
    q=np.exp(-(d**2)/A)
    a1=np.exp(-(a**2)/A)
    r=mpmath.jtheta(3, p*1j, q)
    qr=2*a1*r
    qr=float(str(qr.real))
    return qr/(A**pow)

Best Answer

Both the summation routine I coded to verify the hypothesis and the hypothesis itself was wrong.

Firstly, the summation that I reach from (1) to (2) is correct. However, transition from (2) to (3) is not. $\vartheta_3(.)$ is defined as

$$ \vartheta_3(q,z)=1+2\sum\limits_{i=1}^\infty q^{i^2}\cos(2iz). $$

Rearranging the terms here, we have

$$ \sum\limits_{i=1}^\infty q^{i^2}\cos(2iz)=\frac{\vartheta_3(q,z)-1}{2} $$

The summation limit starts at $i=1$. Thus, (1) becomes

$$ \sum\limits_{i=0}^\infty \frac{1}{A^x}\left[\exp\left(-\frac{(id-a)^2}{A}\right)+\exp\left(-\frac{(id+a)^2}{A}\right)\right]=\frac{2}{A^x}\exp\left(-\frac{a^2}{A}\right)\left[\frac{1}{2}\left(\vartheta_3\left(q,pj\right)-1 \right)+1\right], $$ where $$ q=\exp\left(-\frac{d^2}{A}\right),\phantom{0000} p=\frac{ad}{A}. $$

Note that the last term comes from $i=0$ of the summation.

The verification code can be corrected as follows

def vartheta(A,d,a,pow):
    p=(a*d)/A
    q=np.exp(-(d**2)/A)
    a1=np.exp(-(a**2)/A)
    r=mpmath.jtheta(3, p*1j, q)
    qr=a1*(r+1)
    qr=float(str(qr.real))
    return qr/(A**pow)

As a final note, the first term of the series dominates the summation for small $A$. This is why the result looks perfect for small $A$. As $A$ increases, the contribution of $i=0$ term in the total sum drops, and the result gets more off.

Related Question