Зависящее от времени одномерное уравнение Шредингера с использованием Numpy и SciPy solve_ivp - PullRequest
Новая
       43

Зависящее от времени одномерное уравнение Шредингера с использованием Numpy и SciPy solve_ivp

1 голос
/ 13 апреля 2019

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

enter image description here

* 1007Скажем, у меня есть N пространственных точек (x_i идет от 0 до N-1), и предположим, что мой промежуток времени равен K временным точкам.

Я стремлюсь получить матрицу K by 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:\Users\Pasha\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 704, in runfile
    execfile(filename, namespace)

  File "C:\Users\Pasha\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.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:\Users\Pasha\Anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 454, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)

  File "C:\Users\Pasha\Anaconda3\lib\site-packages\scipy\integrate\_ivp\radau.py", line 288, in __init__
    self.f = self.fun(self.t, self.y)

  File "C:\Users\Pasha\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 139, in fun
    return self.fun_single(t, y)

  File "C:\Users\Pasha\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 21, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)

TypeError: 'numpy.ndarray' object is not callable

В более общем замечании, как я могу заставить python решать Node без определения каждого из них вручную?

Я хочу иметь большой вектор, называемый xdot, где каждая ячейка в этом векторе будет функцией некоторых X [i], и мне кажется, что я не могу этого сделать?или может мой подход совершенно неверный?

Также я чувствую, что, возможно, подключен аргумент "Vectorized" в ivp_solve, но я не понимаю объяснения в документации SciPy.

1 Ответ

1 голос
/ 14 апреля 2019

Проблема, вероятно, в том, что 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)

Этот код решил проблему для меня.

Для краткости вы можете просто написать встроенную функцию, используя лямбду: 

psi= solve_ivp(lambda t,psi : ode_system(psi,t,delta_x,N),[0,15],psi0,method='Radau',max_step=0.1)
...