# Oscillating eigenvectors for 2d-laplace operator

I am trying to calculate the eigenvectors of a Laplace-operator in 2d, with boundary conditions equal to $$u=0$$ if $$x, y$$ are outside of a rectangle defined as $$(0.5, 0.5), (1.5, 1.5)$$. For that I used the following code:

#!/usr/bin/env python3.9

import numpy as np
import scipy as sc
from scipy.sparse import diags, linalg
from scipy import linalg
from scipy import sparse as scsp
from numpy import vectorize
import matplotlib.pyplot as plt
from enum import Enum, IntEnum
import itertools
import sys

np.set_printoptions(threshold=np.inf, linewidth=200)

class scale_dir(IntEnum):
default = 0
up = 1
down = -1

class point_location(IntEnum):
outside = -1
on_line = 0
inside = 1

class point:
def __init__(self, x, y):
self.x = x
self.y = y

return point(self.x + p1.x, self.y + p1.y)

def __sub__(self, p1):
return point(self.x - p1.x, self.y - p1.y)

def __str__(self):
return "(" + str(self.x) + ", " + str(self.y) + ")"

def __equ__(self, p1):
if self.x == p1.x and self.y == p1.y:
return True
else:
return False

def __nequ__(self, p1):
if self.x != p1.x or self.y != p1.y:
return True
else:
return False

class rectangle:
def __init__(self, p0, p1):
self.p0 = p0
self.p1 = p1
if p1.x < p0.x and p0.y < p1.y:
p0.x, p1.x = p1.x, p0.x
if p1.y < p0.y and p0.x < p1.x:
p0.y, p1.y = p1.y, p0.y

def is_point_inside(self, target_point):
if self.p1.x > target_point.x > self.p0.x and self.p1.y > target_point.y > self.p0.y:
return True
else:
return False

def main():
num_grid_points = 64
x_vec = np.linspace(0, 2, num_grid_points)
dx = x_vec - x_vec
y_vec = np.linspace(0, 2, num_grid_points)
dy = y_vec - y_vec

z_matrix = np.zeros((num_grid_points, num_grid_points))

square = rectangle(point(0.5, 0.5), point(1.5, 1.5))

for i in range(num_grid_points):
for j in range(num_grid_points):
z_matrix[i, j] = square.is_point_inside(point(x_vec[i], y_vec[j]))

print(z_matrix)
z_vector = np.reshape(z_matrix, -1)

empty_target_matrix = np.zeros((num_grid_points * num_grid_points, num_grid_points * num_grid_points))

#Test iterating over matrix + matrix-formation
#stencil pattern
for i in range(num_grid_points * num_grid_points):
stencil_vector = np.zeros((2 * num_grid_points + 1))
stencil_vector = 1
stencil_vector[num_grid_points - 1] = 1
stencil_vector[num_grid_points] = -4
stencil_vector[num_grid_points + 1] = 1
stencil_vector[-1] = 1
stencil_vector *= 1 / (dx * dx)
#apply general boundary conditions
insert_vector = stencil_vector.copy()
if i < num_grid_points:
stencil_vector[:num_grid_points - i] = 0
if i > (num_grid_points * num_grid_points - num_grid_points - 1):
stencil_vector[num_grid_points * num_grid_points - i - num_grid_points - 1:] = 0
if i < num_grid_points:
insert_vector = stencil_vector[num_grid_points - i:].copy()
if i > (num_grid_points * num_grid_points - num_grid_points - 1):
insert_vector = stencil_vector[:(num_grid_points * num_grid_points - i - num_grid_points - 1)].copy()
# #Put into matrix
target_row = i
target_col = i - (int)(0.5 * len(stencil_vector))
if target_col < 0:
target_col = 0
for col, elem in enumerate(insert_vector):
empty_target_matrix[i, col + target_col] = elem

#Delete rows and cols, to add boundary conditions
for i, entry in enumerate(z_vector):
if not entry:
empty_target_matrix[i, :] = 0
empty_target_matrix[:, i] = 0
dxdy_matrix = scsp.csr_matrix(empty_target_matrix)

eigenvalues, eigenvectors = scsp.linalg.eigs(dxdy_matrix)

print(eigenvalues)
plt.plot(np.real(eigenvectors[:, 0]))
plt.plot(np.real(eigenvectors[:, 1]))
plt.plot(np.real(eigenvectors[:, 2]))
plt.show()
plt.imshow(np.outer(np.real(eigenvectors[:, 0]), np.real(eigenvectors[:, 0])), extent=[0, 2, 0, 2])
plt.show()

if __name__ == '__main__':
main()


My first three eigenvectors are: i.e. looking similar to the results given in Correct eigenfunctions of Laplace operator by Finite Differences. Nevertheless, I am not able to apply the given solution there (i.e. choosing only the lowest eigenvalues). If I use that approach, I get i.e. something else is going wrong. But I am not able to see what is going wrong, and how I can obtain more correct eigenvectors from the 2d-laplace. Is there an issue with the code, or in the general approach?