Стандартный код решения трехдиагональной матрицы, используемый в al oop, расходится с решениями numpy .linalg.solve () после> 50 итераций.

Question

Я пытаюсь ускорить программу, которую я написал, и я пытался сделать пользовательское решение трехдиагональной матрицы.

Я использую Python 3.8.2 numpy 1.18.2

Я заметил расхождение в моих ответах от правильных, которые оказали значительное влияние на мои результаты примерно за 20 итераций. В попытке исправить отклонение: я вынудил dtype к complex128 и complex64 для массивов / списков (complex64, очевидно, работает хуже). Я превратил все списки в массивы, чтобы сравнить обработку без влияния на результаты. Я пытался без / =, + =, et c. примечание, чтобы избежать каких-либо проблем с выделением памяти, не влияющих на результаты. Я воссоздал методологию декомпрессии LU, используемую linalg.solve (), с сокращенным вычислением трехдиагональности (код не показан) (значительное расхождение произошло при гораздо меньшем количестве итераций)

linalg.solve () обеспечивает правильный ответ для программы, над которой я работаю, и расхождение в моем коде приводит к ошибочным результатам.

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

Вот код, который создает A и b для Ax = b и выполняет мою функцию basi c trislv (), за которой следуют результаты расхождения в одной точке (средней точке) вектора N = 100 после 100 итераций.

Спасибо за взгляд! Мне просто любопытно на этот раз.

import numpy as np
from cmath import exp
import scipy.constants as sc
import matplotlib.pyplot as plt

################################################################
###build tridiagonal solver#################3
def trislv(A,b):

    N=A.shape[1]
    x=[complex(0)]*N #list for solution

    for i in range(N-1):

        A[i,i+1]/=A[i,i]
        b[i]/=A[i,i]
        A[i+1,i+1]-= A[i,i+1]*A[i+1,i]
        b[i+1]-= b[i]*A[i+1,i]
        A[i,i]/=A[i,i]
        A[i+1,i]=0

    b[N-1]/=A[N-1,N-1]
    A[N-1,N-1]/=A[N-1,N-1]

    x[-1]=b[-1]
    for j in range(N-2,-1,-1):
        x[j]=b[j]-A[j,j+1]*x[j+1]
    return x


#################################################################

###Build psi0 and A to test solver##################################

###Physical Parameters##########
m=9.109e-31
L=1e-8
t=0
x0=L/2
s=1e-10 #s for sigma
K=5e10
###Computational Parameters#########
N=1000
a=L/(N+2)
h=1e-18



###Initialize lists for psi(0) and v=B*psi(t) from A*psi(t+h)=B*psi(t)#######################
v=[0.0]*(N+2)
psi0=[0.0]*(N+2)

###building A matrix coefficients###########
a1=1+h*(1j*sc.hbar)/(2*m*a**2)
a2=-h*(1j*sc.hbar)/(4*m*a**2)

###initialize A matrix##########################
A=np.zeros([N,N],complex)

###Build A matrix in row format##############
for i in range(N):
    if i-1>=0:
        A[i,i-1]=a2
    if i+1<=N-1:
        A[i,i+1]=a2
    A[i,i]=a1

###Function to build psi0##############################
def psi(x):
    return exp(-((x-x0)**2)/(2*s**2))*exp(1j*K*x)

###Build psi0 by calling psi()########################
for i in range(1,N+1):
    psi0[i]= psi(a*i)
b=psi0[1:-1]
c=psi0[1:-1]

#############################################################
###test trislv vs. linalg#####################################
for i in range(100):
    y=np.linalg.solve(A,b)
    x=trislv(A,c)

    b=y
    c=x
#x=np.array(x)
print('Linalg.solve solution(midpoint): '+str(np.real(y[500]))+' + '+ str(np.imag(y[500]))+'j')
print(type(y[500]))

print('trislv function solution(midpoint): '+str(np.real(x[500]))+' + '+ str(np.imag(x[500]))+'j')
print(type(x[500]))
d=x[500]-y[500]

print('Difference at midpoint: '+str(np.real(d))+' + '+ str(np.imag(d))+'j')to 

Linalg.solve solution(midpoint): -0.3928294829613219 + -0.42311460439708143j
<class 'numpy.complex128'>
trislv function solution(midpoint): -0.39282994388046255 + -0.42311521934161656j
<class 'numpy.complex128'>
Difference at midpoint: -4.609191406323987e-07 + -6.149445351266714e-07j