0
$\begingroup$

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

    def __add__(self, p1):
        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[1] - x_vec[0]
    y_vec = np.linspace(0, 2, num_grid_points)
    dy = y_vec[1] - y_vec[0]

    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[0] = 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:

enter image description here

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

enter image description here

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?

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.