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$
5
  • 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, 2014 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, 2014 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, 2014 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, 2014 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$ Feb 20, 2014 at 1:54

1 Answer 1

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 agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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