Finite difference discretization on a circle

Question

I am trying to discretize the differential operator $\frac{d^2}{dx^2}$ acting on $S^1 = [0,1]$ using finitely many points around a circle at $0, \frac{1}{N}, \frac{2}{N}, \dots, \frac{N-1}{N}$.

Here is the matrix for $N=5$:

$$\left[ \begin{array}{cc} -2 & 1 & 0 & 0 & 1 \\ 1 & -2 & 1 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 1 & -2 & 1 \\\ 1 &0 & 0 & 1 & -2\end{array} \right]$$

If we solved the equation $\frac{d^2 f}{dx^2} = \lambda f$ we should get $f(x) = e^{\lambda x}$ with $\lambda = 2\pi k$ and $k^2 \in \mathbb{Z}$.

---

When I implemented this in numpy, the reascaling seems all wrong. Can someone explain a little bit of the theory behind this?

N = 100
A = np.zeros((N,N))
d = np.arange(N)
A[d,d] = -2
A[d,(d+1)%N]=1
A[d,(d-1)%N]=1
A

The resulting matrix is like I said. Here $N=100$.

array([[-2.,  1.,  0., ...,  0.,  0.,  1.],
       [ 1., -2.,  1., ...,  0.,  0.,  0.],
       [ 0.,  1., -2., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ..., -2.,  1.,  0.],
       [ 0.,  0.,  0., ...,  1., -2.,  1.],
       [ 1.,  0.,  0., ...,  0.,  1., -2.]])

I came up with this adhoc scaling for the eigenvalues of $-A$ dividing them by $(2 \pi N)^2$.

L = np.linalg.eig(-A)[0]/(2*np.pi)**2*N**2
L = np.round(L)
np.sort(L).astype(int)

The result looks pretty close. I get the sequence of perfect squares. Does that mean $\frac{1}{2\pi N}A$ is the correct rescaling?

array([   0,    1,    1,    4,    4,    9,    9,   16,   16,   25,   25,
         36,   36,   48,   48,   63,   63,   79,   79,   97,   97,  116,
        116,  137,  137,  160,  160,  184,  184,  209,  209,  235,  235,
        263,  263,  291,  291,  320,  320,  350,  350,  381,  381,  412,
        412,  443,  443,  475,  475,  507,  507,  538,  538,  570,  570,
        602,  602,  633,  633,  663,  663,  693,  693,  722,  722,  751,
        751,  778,  778,  804,  804,  830,  830,  853,  853,  876,  876,
        897,  897,  916,  916,  934,  934,  951,  951,  965,  965,  978,
        978,  988,  988,  997,  997, 1004, 1004, 1009, 1009, 1012, 1012,
       1013])

Answer 1

The expressions you are using for eigenvalues and eigenfunctions are wrong (as per Wolfgang Bangerth's comment); therefore the results you are getting are not meaningful at all.

There are analytic expressions for eigenvalues for both continuous and discrete cases for $d^2/dx^2$ on a circular segment – periodic BCs.

For a continuous case on $S^1=[0,L]$, $j$th eigenvalue $\lambda_j$:

$$ \lambda_j^{\text{cont}} = \begin{cases} -\frac{\pi^2}{L^2}j^2,&\quad j \text{ is even}\\ -\frac{\pi^2}{L^2}(j-1)^2,&\quad j \text{ is odd} \end{cases} $$

For a discrete case with $N$ points distributed uniformly around a circle and central differencing scheme:

$$ \lambda_j^{\text{disc}} = \begin{cases} -\frac{4}{h^2}\sin^2\left(\frac{\pi j^2}{2N}\right),&\quad j \text{ is even}\\ -\frac{4}{h^2}\sin^2\left(\frac{\pi (j-1)^2}{2N}\right),&\quad j \text{ is odd} \end{cases} $$

which is effectively the analytical expression for the eigenvalues of your matrix $A$. Here, $h=L/N.$ for a circular segment $S^1$.

I wrote a simple Python script that compares the eigenvalues of the continuous operator and discrete (obtained from the actual matrix using numpy.linalg.eigvals).

For $N=5$, the relative error for the first eigenvalues was:

[0.1248598 0.1248598 0.4272133 0.4272133] (first eigenvalue is 0 and is discarded, eigenvalues have multiplicity 2)

For $N=100$,

[0.00032894 0.00032894 0.00131525 0.00131525]

So, the scheme "converges".

#!/usr/local/bin/python
import numpy as np

# defining problem
L = 1.      # length of the circular interval S = [0,L] (starting is assumed to be fixed
N = 100       # Number of points discretizing circular interval S
h = L/N     # FD step-size (note, it is a circle; thus, N vs. N-1

# forming discretized d^2 / dx^2 operator
A = np.zeros((N, N))
d = np.arange(N)
A[d, d] = -2
A[d, (d+1) % N] = 1
A[d, (d-1) % N] = 1
A *= 1./(h**2)

# calculating eigenvalues
A_eig = np.linalg.eigvals(A)
A_eig = np.flip(np.sort(A_eig), 0)  # decreasing order

# comparing to analytical expressions for eigenvalues in continuous and discrete cases
eig_cont_anal = np.zeros(N)     # continuous problem with periodic BC
eig_disc_anal = np.zeros(N)     # discrete problem with periodic BC

for j in range(1, N+1):
    if j % 2 == 0:  # even eigenvalue
        arg = j
    else:           # odd eigenvalue
        arg = j-1
    eig_cont_anal[j-1] = arg**2
    eig_disc_anal[j-1] = np.sin(np.pi*arg/(2*N))**2

# scaling eigenvalues according to the corresponding analytical expressions
factor_cont = -(np.pi/L)**2
factor_disc = -(4./h**2)
eig_cont_anal *= factor_cont
eig_disc_anal *= factor_disc

# will calculate the relative error between the eigenvalues (except first lambda_0=0) obtained from
# a) analytic eigenvalues for the continuous  problem
# b) numerical eigenvalues for discrete problem with N points

error = np.zeros(N-1)
for j in range(1, N):
    error[j-1] = abs(A_eig[j]-eig_cont_anal[j])/abs(eig_cont_anal[j])

# print(eig_cont_anal); print(A_eig)
K = 4  # number of eigenvalues to print out a relative error
print(error[0:K])