Новичок
Профиль
Группа: Участник
Сообщений: 15
Регистрация: 21.4.2012
Репутация: нет
Всего: нет
|
Здравствуйте, дорогое сообщество Vingrad. Я уже неделю пытаюсь решить одну лабораторную по теме "Одномерное стационарное уравнение Шрёдиннгера в прямоугольной потенциальной яме с бесконечными стенками". Алгоритм Имеем U(x)=V0*v(x) и начальное приближение E=U_min + б, где б>0, U_min - наименьшее значение потенциальной ф-ции U(x). Необходимо найти основное, первое и второе возбуждённые состояния электрона. Делаем это вычислением f через Метод Нумерова (интегрированием вперёд psi(x) и интергрированием назад phi(x) ) и сдвигом E с небольшим шагом. Продолжаем это до тех пор, пока не найдём место. где функция f меняет знак (f1*f2<=0). После этого, при помощи метода бисекций уточняем найденное f2 и строим графики U(x) [жёлтый], phi(x) [синий] , psi(x) [красный], x [зелёный] и указываем узел сшивки(задают произвольно). Суть в том, что программа расчитывает неверно, а где ошибка - понять не могу. Изначальный код выдал преподаватель для встраивания метода бисекции и поиска "перегиба". Дополнительную справку прикладываю в методичке. | Код |
# import time
import math
import numpy as np
# from numpy import *
# from scipy.integrate import odeint
import matplotlib.pyplot as plt
# потенциальная ф-я
def U(x):
return float(V0*(0.5*x*x-1) if abs(x)<L else W)
#return float((-1 + (x*x+L*L)/(2*L*L)) if abs(x) < L else W)
#return float(U0 if abs(x) < L else W)
# ф-я (13)
def q(e, x):
return 2.0*(e-U(x))
# ф-я для вычисления производной по ф-ле (19)
def deriv(Y, h, m):
return (Y[m-2]-Y[m+2]+8.0*(Y[m+1]-Y[m-1]))/(12.0*h)
# метод Нумерова
def num(e, r, n):
# print("r =", r)
F = np.array([c*q(e, X[i]) for i in np.arange(n)])
# print("F =", F)
Psi[0] = 0.0
Fi[n-1] = 0.0
Psi[1] = d1
Fi[n-2] = d2
# интегрирование "вперед" по ф-ле (16)
for i in np.arange(1, n-1, 1):
p1 = 2.0*(1.0 - 5.0*F[i])*Psi[i]
p2 = (1.0 + F[i-1])*Psi[i-1]
Psi[i+1] = (p1 - p2)/(1.0 + F[i+1])
# интегрирование "назад" по ф-ле (17)
for i in np.arange(n-2, 0, -1):
# print("Fi[i+1]=", Fi[i+1], " Fi[i]=", Fi[i])
f1 = 2.0*(1.0 - 5.0*F[i])*Fi[i]
f2 = (1.0 + F[i+1])*Fi[i+1]
Fi[i-1] = (f1 - f2)/(1.0 + F[i-1])
# поиск максимального по модулю эл-та списка Psi
p1 = np.abs(Psi).max()
p2 = np.abs(Psi).min()
big = p1 if p1 > p2 else p2
# масштабирование эл-ов списка Psi
for i in np.arange(0, n-1, 1):
Psi[i] = Psi[i]/big
# масштабирование эл-ов списка Fi так, чтобы F[r]=Psi[r]
coef = Psi[r]/Fi[r]
for i in np.arange(0, n-1, 1):
Fi[i] = coef*Fi[i]
# вычисление f по ф-ле (18) в узле сшивки r
f = deriv(Psi, h, r) - deriv(Fi, h, r)
return f
# границы ямы с а.е.
L = 2.0/0.5292
A = -L
B = +L
# кол-во узлов сетки - нечетное число
n = 5001
# шаг сетки
h = (B-A)/(n-1)
# константа для метода Нумерова
c = h**2/12.0
V0=10
# минимальное значение потенциала а а.е.
U0 = V0*-1.0/27.212
# максимальное значение потенциала на графике
W = V0*2.0/27.212
# описание и инициализация списков
Psi = np.zeros(n)
Fi = np.zeros(n)
F = np.zeros(n)
Psi2 = np.zeros(n)
X = np.linspace(A, B, n)
# задание узла сшивки
r = (n-1)/2+15
d1 = 1.e-9
d2 = d1
# ввод пристрелочного значения энергии с клавиатуры
e = float(input("Начальное приближение E = U_min + "))
e=e+U0
print("e =", e)
f = num(e, r, n)
print("f =", f)
# ввод шага ΔЕ
dE = float(input("ΔE = "))
maxsost = int(input("До какого состояния считать? (0 - основное): "))
maxsost=maxsost+1
f1 = num(e, r, n)
e = e+dE
f2 = num(e, r, n)
step = 0
cursost=0
# метод бисекции
while cursost<maxsost:
while f1*f2>0:
f1 = f2
e = e+dE
f2 = num(e, r, n)
step = step+1
cursost=cursost+1
if cursost<maxsost:
e = e+dE
f2 = num(e, r, n)
exactness = 0.0001
f3 = f2
e1 = e - dE
e2 = e
it=0
e3=e
maxit=50
while abs(f3) > exactness and it<maxit:
e3 = (e2 + e1)/2
f3 = num(e3, r, n)
step = step + 1
if f1*f3<0:
e2 = e3
f2 = f3
else:
e1 = e3
f1 = f3
it += 1
f = f3
print("Расчетное значение e =", e3)
print("f(e) =", f)
i=0
FPl = np.zeros(n)
while i<n:
FPl[i] = Psi[i]*Psi[i]
i = i+1
#print (FPl)
i = 1
h = (B-A)/n
sum = 0
while i<n:
sum = sum + (FPl[i]+FPl[i-1])/2*h
i = i+1
print ("интеграл: ", sum)
i=0
#kor = (math.sqrt(sum))
#print ("kor = ", kor)
while i<n:
FPl[i] = FPl[i] / sum
i = i+1
#print (FPl)
i = 1
sum1 = 0
while i<n:
sum1 = sum1 + (FPl[i]+FPl[i-1])/2*h
i = i+1
print ("после нормировки: ", sum1)
plt.style.use('dark_background')
plt.figure(1, figsize=(12, 8))
plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
plt.rc('text', usetex=True)
plt.rc('font', size=16)
Upot = np.array([U(X[i]) for i in np.arange(n)])
plt.plot(X[1:n-1], Psi[1:n-1], 'r-', linewidth=4.0, label=r"$\psi(x)$")
plt.plot(X[1:n-1], Fi[1:n-1], 'b-', linewidth=2.0, label=r"$\phi(x)$")
plt.plot(X[1:n-1], FPl[1:n-1], 'g-', linewidth=2.0, label=r"$\fpl(x)$")
plt.plot(X, Upot, 'y-', linewidth=6.0, label=r"$U(x)$")
# plt.plot(X, Upot, 'k-', linewidth=6.0, label=r"$U(x)$")
plt.xlabel(r"$x$", fontsize=26, color="w")
plt.ylabel(r"${\psi(x), \phi(x), U(x)}$", fontsize=26, color="w")
plt.grid(True)
plt.legend(fontsize=16, shadow=True, fancybox=True)
plt.plot([X[r]], [Psi[r]], color='green', marker='o', markersize=10)
string1 = "E = " + str(e3)
string2 = "f(E) = " + str(f)
plt.text(-1.5, 2.6, string1, fontsize=18, color='yellow')
plt.text(-1.5, 2.3, string2, fontsize=18, color="yellow")
# we can make labels and lines bigger/thicker by
# making the plot smaller!
F = plt.gcf() # get current figure
size = F.get_size_inches()
# print('size =', size)
fac = 0.9
F.set_size_inches(size[0]*fac, size[1]*fac)
plt.savefig('schrodinger-1.pdf', dpi=300) # save to file
plt.show()
|
MolDiz1a.pdf 442,61 Kb |
[Python] Уравнение Шрёдиннгера в потенциальная яма, Возможная потеря точности 
Загрузка ...
