2
$\begingroup$

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])
$\endgroup$
  • 4
    $\begingroup$ If your matrix is actually circulant, I think you're missing a $1$ in the top right element in the first equation, no? $\endgroup$ – horchler Feb 18 '14 at 18:54
  • $\begingroup$ Indeed. With the last row you are feeding the last element with the first one, but with the scheme you are using you also have to feed the first element with the last one, in a symmetric way. The operator should be symmetric. So the $1$ in the top right as @horchler said. $\endgroup$ – sebas Feb 18 '14 at 19:09
  • 1
    $\begingroup$ And the order of operations in L = np.linalg.eig(-A)[0]/(2*np.pi)**2*N**2 does not match up with dividing by $(2 \pi N)^2$. $\endgroup$ – horchler Feb 18 '14 at 19:11
  • $\begingroup$ I think that you are not dividing by $(2 \pi N)^2$, but dividing by $(2\pi)^2$ and multiplying by $N^2$. You should do L = np.linalg.eig(-A)[0]/(2*np.pi*N)**2 $\endgroup$ – sebas Feb 18 '14 at 23:01
  • 2
    $\begingroup$ Your eigenvalues and eigenfunctions are wrong. First, the function has to be periodic, so you need to have an imaginary exponent, $f(x)=e^{i2\pi k x}$. Second, your eigenvalues are then $\lambda=-(2\pi k)^2$. $\endgroup$ – Wolfgang Bangerth Feb 20 '14 at 1:54
4
$\begingroup$

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])
$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.