[VA] SRC#001 - Spiky Integral

[Thomas Klemm]

(07-13-2018 06:52 PM)ijabbott Wrote:  It's a bit tedious to continue

We can use:
\[cos(nx)=\frac{e^{inx}+e^{-inx}}{2}\]

Thus:
\[ \begin{eqnarray}
cos(mx)cos(nx) &=&(\frac{e^{imx}+e^{-imx}}{2})(\frac{e^{inx}+e^{-inx}}{2}) \\
&=& \frac{e^{i(m+n)x}+e^{-i(m+n)x}}{4}+\frac{e^{i(m-n)x}+e^{-i(m-n)x}}{4} \\
&=&\frac{cos(m+n)}{2}+\frac{cos(m-n)}{2}
\end{eqnarray} \]

Let's forget the factor \(\frac{1}{2}\) for a moment and define:
\[a_{k}=a^{k}+a^{-k}\]

For the same reason as above we have:
\[ \begin{eqnarray}
a_{m}a_{n}&=&(a^{m}+a^{-m})(a^{n}+a^{-n}) \\
&=&a^{m+n}+a^{-(m+n)}+a^{m-n}+a^{-(m-n)} \\
&=&a_{m+n}+a_{m-n}
\end{eqnarray} \]

But of course \(a_{k}=a_{-k}\).

This allows us to rewrite it as:
\[a_{m}a_{n}=a_{m+n}+a_{|m-n|}\]

We want to calculate the product:
\[p_N=\Pi_{k=1}^{N}a_{k}\]

Let's assume we already have \(p_{N-1}\) represented as a sum of \(a_k\) with coefficients \(b_k\):
\[p_{N-1}=\Sigma_{k=1}^{M}b_{k}a_{k}\]

Then \(p_{N}=a_{N}p_{N-1}\) and thus:
\[ \begin{eqnarray}
p_{N} &=& a_{N}\Sigma_{k=1}^{M}b_{k}a_{k} \\
&=& \Sigma_{k=1}^{M}b_{k}a_{N}a_{k} \\
&=& \Sigma_{k=1}^{M}b_{k}(a_{N+k}+a_{N-k})
\end{eqnarray} \]

This Python program allows us to calculate the coefficients \(b_k\):

Code:

def spiky(N):
    b = [0] * N
    b[0] = [0]
    b[1] = [0, 1]
    for k in range(2, N):
        M = len(b[k-1])
        b[k] = [0] * (M + k)
        for i in range(M):
            b[k][k + i] += b[k - 1][i]
            b[k][abs(k - i)] += b[k - 1][i]
    return b

Here's the result for the first couple of values:

Code:

>>> for b in spiky(11):
...   print b
...
[0]
[0, 1]
[0, 1, 0, 1]
[1, 0, 1, 0, 1, 0, 1]
[1, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]
[0, 3, 0, 3, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]
[0, 5, 0, 5, 0, 4, 0, 4, 0, 4, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]
[4, 0, 8, 0, 8, 0, 7, 0, 7, 0, 6, 0, 5, 0, 5, 0, 4, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]
[7, 0, 13, 0, 13, 0, 13, 0, 12, 0, 11, 0, 10, 0, 9, 0, 8, 0, 7, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]
[0, 23, 0, 23, 0, 22, 0, 21, 0, 21, 0, 19, 0, 18, 0, 17, 0, 15, 0, 13, 0, 12, 0, 10, 0, 9, 0, 8, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]
[0, 40, 0, 39, 0, 39, 0, 38, 0, 36, 0, 35, 0, 33, 0, 31, 0, 29, 0, 27, 0, 24, 0, 22, 0, 20, 0, 17, 0, 15, 0, 13, 0, 11, 0, 10, 0, 8, 0, 6, 0, 5, 0, 4, 0, 3, 0, 2, 0, 2, 0, 1, 0, 1, 0, 1]

But we're only interested in the coefficient with the index 0:

Code:

>>> for b in spiky(40):
...   print b[0]
...
0
0
0
1
1
0
0
4
7
0
0
35
62
0
0
361
657
0
0
4110
7636
0
0
49910
93846
0
0
632602
1199892
0
0
8273610
15796439
0
0
110826888
212681976
0
0
1512776590

Now it's time to deal with the factor \(\frac{1}{2}\) that we neglected.
But that's trivial. We just have to add it to each factor \(a_{k}\).
This leads to: \(\frac{1}{2^{N-1}}\)

And since we integrate the constant over \(2\pi\) we loose another 2.

These are the examples for \(N=39\) and \(N=71\):

Code:

>>> b = spiky(100)
>>> print b[39][0], '/', 2**37
1512776590 / 137438953472
>>> print b[71][0], '/', 2**69
2681644149792639400 / 590295810358705651712

However the ratios haven't been reduced.

(07-12-2018 08:10 PM)Gerson W. Barbosa Wrote:  Expand Product{ k=1..n, x^N + 1/x^N } and take the coefficient of the power of x corresponding to the Nth triangular number in the numerator (if there is no correspondence, then the result will be zero). That's your numerator. Your denominator is 2^(N - 1). Multiply the resulting fraction by \(\pi\).

Not sure if I understood that correctly but it might explain the expression:
\[x^N+\frac{1}{x^N}\]

Thanks both for the challenge and the contributions.

Cheers
Thomas