Soluton of an eigenvalue problem using finite differences

eig_prob.py

"""

Problem proposed in:

  https://scicomp.stackexchange.com/a/38968/9667

@author: Nicolás Guarín-Zapata
@date: August 2021
"""

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh, LinearOperator

L = 1
N = 1000
dx = L/(N-1)

A1 = 1.0
A2 = 2.0
A4 = 4.0
A4 = 5.0
D2 = diags([-2., 1., 1.], [0,-1, 1], shape=(N, N))/dx**2
D0 = diags([1.], [0], shape=(N, N))

op1 = LinearOperator((N, N),
                     matvec=lambda v: A1*(D0@v) + A2*(D2@v))


vals, vecs = eigsh(op1, which='SM')
vals_ex = np.array([A1 - A2*k**2*np.pi**2 for k in range(1, 7)])
print("Relative errors:", (-np.sort(-vals) - vals_ex)/vals_ex)