Как применить метод Кранк-Никольсона в питоне к волновому уравнению типа Шредингера - PullRequest
Новая
       66

Как применить метод Кранк-Никольсона в питоне к волновому уравнению типа Шредингера

1 голос
/ 07 июля 2019

Я пытаюсь сделать частицу в имитации коробки без потенциального поля. Мне потребовалось некоторое время, чтобы выяснить, что простые явные и неявные методы нарушают эволюцию унитарного времени, поэтому я прибег к Крэнку-Николсону, который должен быть унитарным. Но когда я пытаюсь это сделать, я обнаруживаю, что это не так Я не уверен, что мне не хватает .. Я использовал следующую формулировку:

formula 1

где T - трехдиагональная матрица Теплица для второй производной по x и

formula 2

Система упрощается до

formula 3

Матрицы A и B:

Я просто решаю эту линейную систему для с помощью разреженного модуля. Математика имеет смысл, и я нашел такую ​​же числовую схему в некоторых работах, что заставило меня поверить, что мой код - вот где проблема.

Вот мой код: 

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import toeplitz
from scipy.sparse.linalg import spsolve
from scipy import sparse

# Spatial discretisation
N = 100
x = np.linspace(0, 1, N)
dx = x[1] - x[0]


# Time discretisation
K = 10000
t = np.linspace(0, 10, K)
dt = t[1] - t[0]

alpha = (1j * dt) / (2 * (dx ** 2))

A = sparse.csc_matrix(toeplitz([1 + 2 * alpha, -alpha, *np.zeros(N-4)]), dtype=np.cfloat)  # 2 less for both boundaries
B = sparse.csc_matrix(toeplitz([1 - 2 * alpha, alpha, *np.zeros(N-4)]), dtype=np.cfloat)

# Initial and boundary conditions (localized gaussian)
psi = np.exp((1j * 50 * x) - (200 * (x - .5) ** 2)) 
b = B.dot(psi[1:-1])
psi[0], psi[-1] = 0, 0

for index, step in enumerate(t):
    # Within the domain
    psi[1:-1] = spsolve(A, b)

    # Enforce boundaries
    # psi[0], psi[N - 1] = 0, 0

    b = B.dot(psi[1:-1])
    # Square integration to show if it's unitary
    print(np.trapz(np.abs(psi) ** 2, dx))

1 Ответ

2 голосов
/ 09 июля 2019

Вы полагаетесь на конструктор Теплица, чтобы получить симметричную матрицу, чтобы записи под диагональю были такими же, как и над диагональю. Тем не менее, документация для scipy.linalg.toeplitz(c, r=None) говорит не "транспонировать" , а

* "Если r не задано, предполагается, что r == сопряжено (c)."

так, чтобы получающаяся матрица была самосопряженной. В этом случае это означает, что записи над диагональю поменялись знаком.

Нет смысла сначала строить плотную матрицу, а затем извлекать разреженное представление. Построить его как разреженную трехдиагональную матрицу с самого начала, используя scipy.sparse.diags

A = sparse.diags([ (N-3)*[-alpha], (N-2)*[1+2*alpha], (N-3)*[-alpha]], [-1,0,1], format="csc");
B = sparse.diags([ (N-3)*[ alpha], (N-2)*[1-2*alpha], (N-3)*[ alpha]], [-1,0,1], format="csc");
...