Persamaan Schrodinger dengan kondisi batas periodik

Saya punya beberapa pertanyaan mengenai hal berikut:

Saya mencoba untuk memecahkan persamaan Schrodinger dalam 1D menggunakan diskritisasi nicolson engkol diikuti dengan membalikkan matriks tridiagonal yang dihasilkan. Masalah saya sekarang telah berkembang menjadi masalah dengan kondisi batas periodik dan karenanya saya telah memodifikasi kode saya untuk menggunakan algoritma Sherman Morrison.

Misalkan vadalah RHS saya pada setiap langkah waktu ketika saya ingin membalikkan matriks tridiagonal. Ukurannya vadalah jumlah titik grid yang saya miliki di atas ruang. Ketika saya mengatur v[0]dan v[-1]dalam hal satu sama lain seperti yang diperlukan dalam situasi periodik saya, persamaan saya meledak. Saya tidak tahu mengapa ini terjadi. Saya menggunakan python2.7 dan scipy's inbuilt resol_banded untuk menyelesaikan persamaan.

Ini mengarahkan saya ke pertanyaan kedua: Saya menggunakan python karena itu adalah bahasa yang saya tahu paling baik, tetapi saya merasa agak lambat (bahkan dengan optimasi yang ditawarkan oleh numpy dan scipy). Saya telah mencoba menggunakan C ++ karena saya cukup terbiasa dengannya. Saya pikir saya akan menggunakan GSL yang akan dioptimalkan BLAS, tetapi tidak menemukan dokumentasi untuk membuat vektor yang kompleks atau memecahkan matriks tridiagonal dengan vektor-vektor bernilai kompleks seperti itu.

Saya ingin objek dalam program saya karena saya merasa itu akan menjadi cara termudah bagi saya untuk menggeneralisasi nanti untuk menyertakan kopling antara fungsi gelombang sehingga saya tetap berpegang pada bahasa berorientasi objek.

Saya bisa mencoba menulis pemecah matriks tridiagonal dengan tangan, tetapi saya mengalami masalah ketika saya melakukannya dengan python. Ketika saya berevolusi berkali-kali dengan langkah waktu yang lebih halus dan lebih halus, kesalahan menumpuk dan memberi saya omong kosong. Dengan mengingat hal ini, saya memutuskan untuk menggunakan metode bawaan.

Setiap saran sangat dihargai.

EDIT: Berikut ini cuplikan kode yang relevan. Notasi tersebut dipinjam dari halaman Wikipedia pada persamaan tridiagonal matrix (TDM). v adalah RHS dari algoritma crank nicolson pada setiap langkah waktu. Vektor a, b dan c adalah diagonal TDM. Algoritma yang dikoreksi untuk case periodik adalah dari CFD Wiki . Saya telah melakukan sedikit penggantian nama. Apa yang mereka sebut u, v Saya telah memanggil U, V (dikapitalisasi). Saya telah memanggil q pelengkap, y solusi sementara dan solusi aktual self.currentState. Penugasan v [0] dan v [-1] adalah apa yang menyebabkan masalah di sini dan karenanya telah dikomentari. Anda dapat mengabaikan faktor-faktor gamma. Mereka adalah faktor non-linear yang digunakan untuk memodelkan Kondensat Bose Einstein.

for T in np.arange(self.timeArraySize):
        for i in np.arange(0,self.spaceArraySize-1):
            v[i] = Y*self.currentState[i+1] + (1-2*Y)*self.currentState[i] + Y*self.currentState[i-1] - 1j*0.5*self.timeStep*potential[i]*self.currentState[i] - self.gamma*1j*0.5*self.timeStep*(abs(self.currentState[i])**2)*self.currentState[i]
            b[i] = 1+2*Y + 1j*0.5*self.timeStep*potential[i] + self.gamma*self.timeStep*1j*0.5*(abs(self.currentState[i])**2)

        #v[0] = Y*self.currentState[1] + (1-2*Y)*self.currentState[0] + Y*self.currentState[-1] - 1j*0.5*self.timeStep*potential[0]*self.currentState[0]# - self.gamma*1j*0.5*self.timeStep*(abs(self.currentState[0])**2)*self.currentState[0]
        #v[-1] = Y*self.currentState[0] + (1-2*Y)*self.currentState[-1] + Y*self.currentState[-2] - 1j*0.5*self.timeStep*potential[-1]*self.currentState[-1]# - self.gamma*1j*0.5*self.timeStep*(abs(self.currentState[-1])**2)*self.currentState[-1]
        b[0] = 1+2*Y + 1j*0.5*self.timeStep*potential[0] + self.gamma*self.timeStep*1j*0.5*(abs(self.currentState[0])**2)
        b[-1] = 1+2*Y + 1j*0.5*self.timeStep*potential[-1] + self.gamma*self.timeStep*1j*0.5*(abs(self.currentState[-1])**2)

        diagCorrection[0], diagCorrection[-1] = - b[0], - c[-1]*a[0]/b[0]

        tridiag = np.matrix([
            c,
            b - diagCorrection,
            a,
        ])

        temp = solve_banded((1,1), tridiag, v)

        U = np.zeros(self.spaceArraySize, dtype=np.complex64)
        U[0], U[-1] = -b[0], c[-1]

        V = np.zeros(self.spaceArraySize, dtype=np.complex64)
        V[0], V[-1] = 1, -a[0]/b[0]

        complement = solve_banded((1,1), tridiag, U)

        num = np.dot(V, temp)
        den = 1 + np.dot(V, complement)

        self.currentState = temp  - (num/den)*complement

Pertanyaan kedua

Kode yang memanggil Scipy / Numpy biasanya hanya cepat jika dapat di-vektor; Anda seharusnya tidak memiliki sesuatu yang "lambat" di dalam lingkaran python. Bahkan kemudian, itu cukup banyak tidak terhindarkan bahwa itu akan setidaknya sedikit lebih lambat daripada sesuatu yang menggunakan perpustakaan serupa dalam bahasa yang dikompilasi.

for i in np.arange(0,self.spaceArraySize-1):
            v[i] = Y*self.currentState[i+1] + (1-2*Y)*self.currentState[i]   ...
            b[i] = 1+2*Y + 1j*0.5*self.timeStep*potential[i] + ...

Ini yang saya maksud dengan "slow in a python loop". Python forsangat lambat untuk sebagian besar aplikasi numerik, dan Scipy / Numpy tidak memengaruhi ini sama sekali. Jika Anda akan menggunakan python, loop dalam ini harus dinyatakan sebagai satu atau dua fungsi Numpy / Scipy, yang mungkin atau mungkin tidak disediakan oleh pustaka tersebut. Jika mereka tidak menyediakan sesuatu yang memungkinkan Anda untuk mengulangi array seperti ini dan mengakses elemen yang berdekatan, python adalah alat yang salah untuk apa yang ingin Anda lakukan.

Juga, Anda melakukan inversi daripada penyelesaian matriks-vektor. Sebuah inversi matriks, diikuti oleh matriks-vektor multiply, jauh lebih lambat daripada penyelesaian matriks-vektor. Ini hampir pasti hal yang memperlambat kode Anda lebih dari apa pun.

Jika Anda ingin menggunakan C / C ++, GSL agak kurang dalam hal aljabar linier kompleks. Saya akan merekomendasikan menggunakan BLAS atau LAPACK secara langsung, atau menggunakan perpustakaan seperti PETSc atau Trilinos. Jika Anda telah menginstal MKL, Anda dapat menggunakannya juga. Anda mungkin juga ingin memeriksa Fortran 2008, yang berorientasi objek.

Matriks Anda jarang, jadi pastikan Anda menggunakan perpustakaan yang jarang.

Saya juga akan mengatakan bahwa apa yang Anda lakukan di sini tampaknya cukup rendah sehingga orientasi objek mungkin tidak menjadi perhatian utama Anda. Fortran 90+ array mungkin cocok dengan yang Anda butuhkan, dan kompiler F90 dapat secara otomatis melakukan paralelisasi beberapa loop.

Juga, Anda mungkin ingin memeriksa Oktaf atau Matlab, yang memiliki sparse()fungsi. Jika digunakan dengan benar, ini harus dapat berjalan cukup cepat.