حل معادلهی شرودینگر ذرهی آزاد#
ویدئوی جلسه
تبدیل فوریه یکی از ابزارهای قدرتمند در ریاضیات و فیزیک است که به ما امکان میدهد توابع را از حوزهی زمان (یا مکان) به حوزهی فرکانس (یا تکانه) تبدیل کنیم. این تبدیل در حل معادلات دیفرانسیل، پردازش سیگنالها و بسیاری از زمینههای دیگر کاربردهای گستردهای دارد. در این بخش، به بررسی تاریخچه، مفاهیم پایه و کاربرد تبدیل فوریه در حل معادلات دیفرانسیل میپردازیم.
تبدیل فوریه توسط ژان-باپتیست ژوزف فوریه، ریاضیدان فرانسوی، در قرن نوزدهم معرفی شد. فوریه در مطالعهی معادلات دیفرانسیل جزئی مربوط به انتقال گرما، متوجه شد که میتوان هر تابع تناوبی را به صورت مجموعی از توابع سینوسی و کسینوسی بیان کرد. این ایده پایهی تبدیل فوریه را تشکیل داد.
مقدمه#
۲. تبدیل فوریه چیست؟#
تبدیل فوریه یک تابع \( f(x) \) را از حوزهی زمان (یا مکان) به حوزهی فرکانس تبدیل میکند. تبدیل فوریهی پیوسته به صورت زیر تعریف میشود:
که در آن:
\( f(x) \): تابع در حوزهی زمان (یا مکان).
\( \hat{f}(\xi) \): تابع در حوزهی فرکانس.
\( \xi \): فرکانس.
تبدیل فوریهی گسسته (DFT) و تبدیل فوریهی سریع (FFT) نیز نسخههای گسستهی این تبدیل هستند که برای دادههای گسسته استفاده میشوند.
۳. کاربرد تبدیل فوریه در حل معادلات دیفرانسیل#
تبدیل فوریه به ویژه در حل معادلات دیفرانسیل خطی با ضرایب ثابت بسیار مفید است. این روش به ما امکان میدهد معادلات دیفرانسیل را در حوزهی فرکانس حل کنیم و سپس با استفاده از تبدیل معکوس فوریه، جواب را به حوزهی زمان بازگردانیم.
الف) مراحل حل معادلات دیفرانسیل با تبدیل فوریه#
تبدیل معادله به حوزهی فرکانس: با اعمال تبدیل فوریه به معادلهی دیفرانسیل، معادلهای در حوزهی فرکانس به دست میآید.
حل معادله در حوزهی فرکانس: معادلهی تبدیلشده معمولاً سادهتر است و میتوان آن را حل کرد.
تبدیل معکوس به حوزهی زمان: با استفاده از تبدیل معکوس فوریه، جواب را به حوزهی زمان بازمیگردانیم.
ب) مثال: حل معادلهی دیفرانسیل موج#
معادلهی دیفرانسیل موج به صورت زیر است:
با اعمال تبدیل فوریه به این معادله، داریم:
این معادله در حوزهی فرکانس سادهتر است و میتوان آن را حل کرد. سپس با تبدیل معکوس فوریه، جواب در حوزهی زمان به دست میآید.
۴. تبدیل فوریهی سریع (FFT)#
تبدیل فوریهی سریع (FFT) یک الگوریتم بهینهشده برای محاسبهی تبدیل فوریهی گسسته (DFT) است. این الگوریتم زمان محاسبات را از \( O(N^2) \) به \( O(N \log N) \) کاهش میدهد و در حل معادلات دیفرانسیل به صورت عددی بسیار مفید است.
کاربردهای FFT#
پردازش سیگنال : تجزیه سیگنالهای صوتی و تصویری به اجزای فرکانسی برای فیلتر کردن، فشردهسازی و تجزیه و تحلیل.
پردازش تصویر: فیلتر کردن، فشردهسازی و تجزیه و تحلیل تصاویر با استفاده از تبدیل فوریه.
ارتباطات: مدولاسیون و دمدولاسیون سیگنالها در سیستمهای ارتباطی.
پزشکی: تجزیه و تحلیل سیگنالهای EEG و ECG برای تشخیص بیماریها.
نجوم: تجزیه و تحلیل سیگنالهای رادیویی از فضا.
حل معادلات دیفرانسیل جزئی: FFT برای حل معادلات دیفرانسیل جزئی در حوزهی فرکانس استفاده میشود.
محاسبهی مشتقات: با استفاده از FFT، میتوان مشتقات توابع را به صورت عددی محاسبه کرد.
امکانات موجود در numpy.fft#
تبدیل فوریه#
توابع زیر یک آرایه را گرفته و تبدیل فوریهی آن را برمیگردانند. دقت کنید که اگر تابع حقیق باشد* باید از rfft استفاده شود:
numpy.fft.fftnumpy.fft.rfft
اگر تابع چند متغیره داشته باشیم، این توابع میتوانند تبدیل فوریه نسبت به هر متغیر دلخواه را حساب کنند. برای توضیحات این امکان به مستندات این کتابخانه مراجعه کنید.
تبدیل معکوس فوریه#
توابع زیر یک آرایه را گرفته و تبدیل فوریهی معکوس آن را برمیگردانند. دقت کنید که اگر تابع حقیق باشد* باید از rfft استفاده شود:
numpy.fft.ifftnumpy.fft.irfft
تبدیل فوریهی چند بعدی#
اگر بخواهید تبدیل فوریهی یک تابع چند متغیره را همزمان نسبت به چند متغیر آن حساب کنید، میتوانید از دستورات زیر استفاده کنید:
numpy.fft.fft2,numpy.fft.ifft2,numpy.fft.rfft2,numpy.fft.irfft2: برای توابع دو بعدیnumpy.fft.fftn,numpy.fft.ifftn,numpy.fft.rfftn,numpy.fft.irfftn: برای بیش از دو بعد
محاسبهی بازهی فرکانسها#
با دانستن بازه ی زمانی و خانه بندی در زمان، میتوان بازهی فرکانسها را به دست آورد. برای این کار توابع زیر را میتوان استفاده کرد:
numpy.fft.fftfreqnumpy.fft.rfftfreq
از آنجا که خروجی این توابع مرتب شده نیستند و اول فرکانسهای مثبت و سپس فرکانسهای منفی قرار دارند، میتوان از تابع fftshit برای مرتب سازی استفاده کرد.
معادلهی شرودینگر ذرهی آزاد#
در مکانیک کوانتومی، معادله شرودینگر یک معادله دیفرانسیلی است که رفتار یک سیستم کوانتومی را توصیف میکند. برای یک ذره آزاد که تحت هیچ نیروی خارجی قرار ندارد، معادله شرودینگر به شکل زیر است:
که در آن:
\( \psi(x,t) \) تابع موج ذره است که احتمال یافتن ذره در موقعیت \( x \) و زمان \( t \) را نشان میدهد.
\( \hbar \) ثابت پلانک تقسیم بر \( 2\pi \) است.
\( m \) جرم ذره است.
حل معادله شرودینگر با استفاده از تبدیل فوریه#
برای حل معادله شرودینگر برای ذره آزاد، از تبدیل فوریه استفاده میکنیم. ابتدا تابع موج \( \psi(x,t) \) را به فضای فرکانس تبدیل میکنیم. تبدیل فوریه تابع موج به شکل زیر است:
که در آن \( \tilde{\psi}(k,t) \) تابع موج در فضای موجی (فضای \( k \)) است و \( k \) عدد موج است.
معادله شرودینگر در فضای \( k \)#
با استفاده از تبدیل فوریه، معادله شرودینگر به فضای \( k \) تبدیل میشود. با مشتقگیری از معادله شرودینگر و اعمال تبدیل فوریه، معادلهای برای \( \tilde{\psi}(k,t) \) به دست میآید:
این معادله یک معادله دیفرانسیلی ساده برای \( \tilde{\psi}(k,t) \) است که میتوان آن را حل کرد. حل عمومی این معادله به شکل زیر است:
که \( \tilde{\psi}(k,0) \) تابع موج اولیه در فضای \( k \) است.
بازگشت به فضای \( x \)#
برای بازگرداندن تابع موج به فضای \( x \)، از تبدیل فوریه معکوس استفاده میکنیم:
رابطه تابع موج در لحظه \( t + \Delta t \)#
برای به دست آوردن تابع موج در زمان \( t + \Delta t \)، میتوانیم از حل به دست آمده برای \( \tilde{\psi}(k,t) \) استفاده کنیم. از آنجایی که تابع موج در فضای \( k \) به صورت نمایی تغییر میکند، داریم:
سپس با استفاده از تبدیل فوریه معکوس، تابع موج در لحظه \( t + \Delta t \) را به دست میآوریم:
این معادله تابع موج را در لحظه \( t + \Delta t \) به دست میدهد.
پیاده سازی#
حالت اولیه#
حالت اولیهی سیستم را میتوانیم با یک تابع گاوسی توصیف کنیم.
import numpy as np
sigma = 10 # set Δx value (σ = 0.1)
p0 = 20 # set initial momentum value
x0 = 0 # set initial position value (x0 = 0)
N = 20000 # set number of points (N = 2000)
L = 100*sigma # set length of the box (L = 100σ)
hbar = 1 # set Planck's constant (ℏ = 1)
x = np.linspace(-L/2 , L/2 , N) # create x array
psi0 = np.exp( -(x-x0)**2 / (4*sigma**2) ) * np.exp(p0*x*1j/hbar) / np.sqrt( np.sqrt(2*np.pi)*sigma)
دقت کنید که تابع بالا پیاده سازی شده ی تابع گاوسی زیر است و مقادیر آن اعداد مختلط هستند:
خواص این تابع صورت است:
\( \int_{-\infty}^{\infty} |\psi_0(x)|^2 dx = 1 \)
\( \langle x \rangle = \int_{-\infty}^{\infty} x|\psi_0(x)|^2 dx = x_0 = 0 \)
\( \langle p \rangle = \int_{-\infty}^{\infty} \psi_0^*(x) \left( -i\hbar \frac{d}{dx} \right) \psi_0(x) dx = p_0 \)
\( \sigma_x = \sqrt{\langle x^2 \rangle - \langle x \rangle^2} = \sigma \)
\( \sigma_p = \sqrt{\langle p^2 \rangle - \langle p \rangle^2} = \frac{\hbar}{2\sigma} \)
\( \sigma_x \sigma_p = \frac{\hbar}{2} \)
میتوانیم این خواص را به راحتی به صورت عددی محاسبه کنیم:
# calculate the normalization of the psi0
dx = x[1] - x[0]
integral = np.sum(np.abs(psi0)**2) * dx
print(f'integral of |psi0|^2 ={integral:.3f}')
integral of |psi0|^2 =1.000
#average and standard deviation of position
avg_x = np.sum(np.abs(psi0)**2 * x) * dx
std_x = np.sqrt( np.sum(np.abs(psi0)**2 * (x-avg_x)**2) * dx )
print(f"average of x ={avg_x:.3f}")
print(f"standard deviation of x ={std_x:.3f}")
average of x =-0.000
standard deviation of x =10.000
# average and standard deviation of momentum
avg_p = np.sum(np.conj(psi0) * (-1j*hbar*np.gradient(psi0,x)) ) * dx
std_p = np.sqrt( np.sum(np.conj(psi0) * (-1*hbar*hbar*np.gradient(np.gradient(psi0,x) , x) )) * dx - avg_p**2 )
print(f"average of p ={np.real(avg_p):.2f}")
print(f"standard deviation of p ={np.real(std_p):.2f}")
print(f"undertainty principle σxσp={np.real(std_p*std_x):.2f}")
average of p =16.83
standard deviation of p =0.03
undertainty principle σxσp=0.27
و همچنینی میتوانیم شکل این تابع را رسم کنیم:
# plot the wave function
import matplotlib.pyplot as plt
range_of_interest_x = np.abs(x-x0)< 4*sigma # range of interest in x space, where the wave function is significant : 4σ around x0
plt.plot(x[ range_of_interest_x ], np.abs(psi0[range_of_interest_x ])**2)
plt.xlabel('x')
plt.ylabel('$|\psi_0|^2$')
plt.title('Initial wave function')
plt.show()
تبدیل فوریه حالت اولیه#
حال میتوانیم تبدیل فوریه ی حالت اولیه را حساب کنیم و با استفاده از آن نیز خواص مومنتومی تابع موج را بررسی کنیم
psi0_k = np.fft.fftshift(np.fft.fft(psi0)) * dx / np.sqrt(2*np.pi) # Fourier transform of psi0
k = np.fft.fftshift(np.fft.fftfreq(N)) *2*np.pi / (dx)
dk = k[1] - k[0]
sigma_p = hbar/(2*sigma)
range_of_interest_p = np.abs(k-p0)<4*sigma_p # range of interest in momentum space, a window of 4σ_p around p0
plt.plot(k[range_of_interest_p], np.abs(psi0_k[range_of_interest_p])**2 )
plt.xlabel('k')
plt.ylabel('$|\psi_0(k)|^2$')
plt.title('Initial wave function in momentum space')
plt.show()
محاسبات مربوط به تکانه نیز در فضای فوریه ساده است:
\( \langle p \rangle = \int_{-\infty}^{\infty} p|\tilde{\psi_0}(p)|^2 dp = p_0 \)
\( \langle p^2 \rangle = \int_{-\infty}^{\infty} p^2|\tilde{\psi_0}(p)|^2 dp = p_0 \)
# calculate the average and standard deviation of momentum in momentum space
size_k = np.sum(np.abs(psi0_k)**2) * dk
avg_k = np.sum(np.abs(psi0_k)**2 * k) * dk
std_k = np.sqrt( np.sum(np.abs(psi0_k)**2 * (k-avg_k)**2) * dk )
print(f"size of psi0_k ={size_k:.3f}")
print(f"average of k ={avg_k:.3f}")
print(f"standard deviation of k ={std_k:.3f}")
size of psi0_k =1.000
average of k =20.000
standard deviation of k =0.050
محاسبه ی تحول زمانی#
برای محاسبه ی مقدار \(\tilde{\psi(k)}\) در زمان \(t+\Delta t\) میتوانیم از رابطهی زیر استفاده کنیم:
میتوانیم این تحول زمانی را به صورت زیر پیاده سازی کنیم:
import tqdm
# implement the time evolution in momentum space
m = 0.001 # set mass (m = 1)
dt = 0.0005 # set time step (dt = 1)
t = 0 # set initial time (t = 0)
T = 5 # set total time (T = 1000)
Nt = int(T/dt) # calculate the number of time steps
times_to_plot = []
states_to_plot_x = []
states_to_plot_k = []
x_avgs = []
x_stds = []
normalizations = []
p_avgs = []
p_stds = []
psi_k = psi0_k # set the initial wave function in momentum space
psi_k /= np.sqrt(size_k) # normalize the wave function in momentum space
times_to_plot.append(0)
states_to_plot_x.append(psi0)
states_to_plot_k.append(psi0_k)
x_avgs.append( np.sum(np.abs(psi0)**2 * x) * dx )
x_stds.append( np.sqrt( np.sum(np.abs(psi0)**2 * (x-x_avgs[-1])**2) * dx ) )
for i in tqdm.tqdm(range(Nt)):
psi_k *= np.exp(-1j * hbar * k*k * dt / (2 * m)) # Apply kinetic energy operator correctly
#psi_k -= psi_k * ( 1j * hbar * k*k * dt / (2*m) ) # apply the kinetic energy operator
#psi_k = psi_k / np.sqrt(size_k) # normalize the wave function in momentum space
p_avgs.append( np.sum(np.abs(psi_k)**2 * k) * dk )
p_stds.append( np.sqrt( np.sum(np.abs(psi_k)**2 * (k-p_avgs[-1])**2) * dk ) )
if i % 200 == 0: # plot the wave function every 100 time steps
psi = np.fft.ifft( np.fft.ifftshift( psi_k) ) * dk / np.sqrt(2*np.pi) # inverse Fourier transform of psi_k
psi /= np.sqrt( np.sum(np.abs(psi)**2) * dx )
times_to_plot.append(t)
states_to_plot_x.append(psi)
states_to_plot_k.append(psi_k)
x_avgs.append( np.sum(np.abs(psi)**2 * x) * dx )
x_stds.append( np.sqrt( np.sum(np.abs(psi)**2 * (x-x_avgs[-1])**2) * dx ) )
t += dt # update the time
0%| | 0/10000 [00:00<?, ?it/s]
100%|██████████| 10000/10000 [00:16<00:00, 603.68it/s]
import matplotlib.animation as animation
from IPython.display import HTML
ims = []
fig = plt.figure()
for state in tqdm.tqdm(states_to_plot_x):
im = plt.plot(x , np.abs( state )**2 , animated=True)
ims.append( im )
ani = animation.ArtistAnimation(fig, ims, interval=500, blit=True,
repeat_delay=1000)
output = HTML(ani.to_jshtml())
plt.close()
output
100%|██████████| 51/51 [00:00<00:00, 251.81it/s]
plt.plot(list(range(len(p_stds))), p_stds , label='standard deviation of p')
plt.plot(list(range(Nt)), p_avgs , label='average of p')
plt.xlabel('time')
plt.legend()
plt.show()
plt.plot(times_to_plot, x_stds , label='standard deviation of x')
plt.plot(times_to_plot, x_avgs , label='average of x')
plt.xlabel('time')
plt.legend()
plt.show()
تکلیف#
تکلیف تبدیل فوریه
در این مسئله، از شما خواسته میشود تا یک انتگرال را با استفاده از تبدیل فوریه سریع (FFT) به دو روش مختلف محاسبه کنید:
انتگرال حاصلضرب تبدیل فوریهی دو تابع.
استفاده از قضیهی کانولوشن.
۱. تعریف مسئله
مسئله: محاسبهی انتگرال وابسته به پارامتر \( a \) با استفاده از تبدیل فوریه سریع (FFT)
در این مسئله، از شما خواسته میشود تا یک انتگرال وابسته به پارامتر \( a \) را با استفاده از تبدیل فوریه سریع (FFT) به دو روش مختلف محاسبه کنید:
انتگرال حاصلضرب تبدیل فوریهی دو تابع.
استفاده از قضیهی کانولوشن.
سپس نتایج را بر حسب پارامتر \( a \) رسم کرده و مقایسه کنید.
۱. تعریف مسئله
فرض کنید دو تابع \( f(x) \) و \( g(x) \) به صورت زیر تعریف شدهاند:
میخواهیم انتگرال زیر را به عنوان تابعی از پارامتر \( a \) محاسبه کنیم:
۲. روش اول: انتگرال حاصلضرب تبدیل فوریهی دو تابع
الف) تبدیل فوریهی توابع
ابتدا تبدیل فوریهی توابع \( f(x - a) \) و \( g(x) \) را محاسبه کنید:
ب) انتگرال حاصلضرب تبدیل فوریهها
با استفاده از قضیهی پارسوال (Parseval’s Theorem)، انتگرال \( I(a) \) را میتوان به صورت زیر محاسبه کرد:
که در آن \( \hat{g}^*(\xi) \) مزدوج مختلط \( \hat{g}(\xi) \) است.
ج) محاسبهی عددی با FFT
توابع \( f(x - a) \) و \( g(x) \) را در یک بازهی مناسب نمونهبرداری کنید.
با استفاده از FFT، تبدیل فوریهی این توابع را محاسبه کنید.
حاصلضرب \( \hat{f}_a(\xi) \cdot \hat{g}^*(\xi) \) را محاسبه و انتگرال آن را با استفاده از روشهای عددی (مانند ذوزنقهای) محاسبه کنید.
۳. روش دوم: استفاده از قضیهی کانولوشن
الف) قضیهی کانولوشن
قضیهی کانولوشن بیان میکند که تبدیل فوریهی کانولوشن دو تابع برابر با حاصلضرب تبدیل فوریهی آنها است:
که در آن \( f * g \) کانولوشن دو تابع \( f \) و \( g \) است.
ب) محاسبهی انتگرال با کانولوشن
با استفاده از قضیهی کانولوشن، انتگرال \( I(a) \) را میتوان به صورت زیر محاسبه کرد:
که در آن \( (f * g)(a) \) مقدار کانولوشن دو تابع در نقطهی \( x = a \) است.
ج) محاسبهی عددی با FFT
توابع \( f(x) \) و \( g(x) \) را در یک بازهی مناسب نمونهبرداری کنید.
با استفاده از FFT، تبدیل فوریهی این توابع را محاسبه کنید.
حاصلضرب \( \hat{f}(\xi) \cdot \hat{g}(\xi) \) را محاسبه و تبدیل معکوس FFT را اعمال کنید تا کانولوشن \( f * g \) به دست آید.
مقدار \( (f * g)(a) \) را به عنوان جواب انتگرال \( I(a) \) در نظر بگیرید.
۴. رسم نتایج و مقایسه
الف) رسم نتایج
مقادیر \( I(a) \) را برای بازهای از مقادیر \( a \) (مثلاً \( a \in [-5, 5] \)) محاسبه کنید.
نتایج حاصل از دو روش را بر حسب \( a \) رسم کنید.
ب) مقایسهی نتایج
دو نمودار را با هم مقایسه کنید و بررسی کنید که آیا نتایج دو روش با هم همخوانی دارند.
اگر اختلافی وجود دارد، دلیل آن را تحلیل کنید (مثلاً خطاهای عددی یا محدودیتهای نمونهبرداری).