时间依赖的1D Schroedinger方程使用Numpy和SciPy solve_ivp
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了时间依赖的1D Schroedinger方程使用Numpy和SciPy solve_ivp相关的知识,希望对你有一定的参考价值。
我试图用有限差分方法求解一维时间依赖的Schroedinger方程,下面是方程看起来如何以及如何进行离散化
假设我有N个空间点(x_i从0变为N-1),并假设我的时间跨度是K个时间点。
我努力获得K×N矩阵。每行(j)将是时间t_j的函数
我怀疑我的问题是我以错误的方式定义了耦合方程的系统。
我的边界条件是psi = 0,或者在框的两侧有一些常数,所以我将X跨度两侧的颂歌设为零
我的代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#Defining the length and the resolution of our x vector
length = 2*np.pi
delta_x = .01
# create a vector of X values, and the number of X values
def create_x_vector(length, delta_x):
x = np.arange(-length, length, delta_x)
N = len(x)
return x, N
# create initial condition vector
def create_initial_cond(x,x0,Gausswidth):
psi0 = np.exp((-(x-x0)**2)/Gausswidth)
return psi0
#create the system of ODEs
def ode_system(psi,t,delta_x,N):
psi_t = np.zeros(N)
psi_t[0]=0
psi_t[N-1]=0
for i in range(1,N-1):
psi_t[i] = (psi[i+1]-2*psi[i]+psi[i-1])/(delta_x)**2
return psi_t
#Create the actual time, x and initial condition vectors using the functions
t = np.linspace(0,15,5000)
x,N= create_x_vector(length,delta_x)
psi0 = create_initial_cond(x,0,1)
psi = np.zeros(N)
psi= solve_ivp(ode_system(psi,t,delta_x,N),[0,15],psi0,method='Radau',max_step=0.1)
运行后我得到一个错误:
runfile('D:/Studies/Project/Simulation Test/Test2.py', wdir='D:/Studies/Project/Simulation Test')
Traceback (most recent call last):
File "<ipython-input-16-bff0a1fd9937>", line 1, in <module>
runfile('D:/Studies/Project/Simulation Test/Test2.py', wdir='D:/Studies/Project/Simulation Test')
File "C:UsersPashaAnaconda3libsite-packagesspyder_kernelscustomizespydercustomize.py", line 704, in runfile
execfile(filename, namespace)
File "C:UsersPashaAnaconda3libsite-packagesspyder_kernelscustomizespydercustomize.py", line 108, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)
File "D:/Studies/Project/Simulation Test/Test2.py", line 35, in <module>
psi= solve_ivp(ode_system(psi,t,delta_x,N),[0,15],psi0,method='Radau',max_step=0.1)
File "C:UsersPashaAnaconda3libsite-packagesscipyintegrate\_ivpivp.py", line 454, in solve_ivp
solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
File "C:UsersPashaAnaconda3libsite-packagesscipyintegrate\_ivp
adau.py", line 288, in __init__
self.f = self.fun(self.t, self.y)
File "C:UsersPashaAnaconda3libsite-packagesscipyintegrate\_ivpase.py", line 139, in fun
return self.fun_single(t, y)
File "C:UsersPashaAnaconda3libsite-packagesscipyintegrate\_ivpase.py", line 21, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
TypeError: 'numpy.ndarray' object is not callable
更一般地说,如何在不手动定义每个和其中一个的情况下使python解决N ode?
我想要一个名为xdot的大向量,这个向量中的每个单元格都是X [i]的函数,我似乎没有做到这一点?或者我的方法完全错了?
我也觉得ivp_solve的“Vectorized”参数可以连接,但我不理解SciPy文档中的解释。
问题可能是solve_ivp
期望第一个参数的函数,并且你提供了ode_system(psi,t,delta_x,N)
,这导致了一个矩阵(因此你得到type error - ndarray
)。
你需要为solve_ivp
提供一个接受两个变量的函数,t
和y
(在你的情况下是psi)。它可以这样做:
def temp_function(t, psi):
return ode_system(psi,t,delta_x,N)
然后,你的最后一行应该是:
psi= solve_ivp(temp_function,[0,15],psi0,method='Radau',max_step=0.1)
这段代码为我解决了这个问题。
对于这样做的简写方法,您也可以使用Lambda编写内联函数:
psi= solve_ivp(lambda t,psi : ode_system(psi,t,delta_x,N),[0,15],psi0,method='Radau',max_step=0.1)
以上是关于时间依赖的1D Schroedinger方程使用Numpy和SciPy solve_ivp的主要内容,如果未能解决你的问题,请参考以下文章