- 相關(guān)推薦
定態(tài)薛定諤方程的MATLAB求解(一)
利用矩陣法對(duì)定態(tài)薛定諤方程的MATLAB求解
摘要：本文首先對(duì)薛定諤方程的提出及發(fā)展做了一個(gè)簡單介紹。然后，以在一維空間運(yùn)動(dòng)的粒子構(gòu)成的諧振子的體系為例，詳細(xì)介紹了矩陣法求解薛定諤方程的過程及公式推導(dǎo)。最后，通過MATLAB編程仿真實(shí)現(xiàn)了求解結(jié)果。
關(guān)鍵詞：定態(tài)薛定諤方程求解 矩陣法 MATLAB仿真
薛定諤方程簡介
1.1背景資料
薛定諤方程是由奧地利物理學(xué)家薛定諤提出的量子力學(xué)中的一個(gè)基本方程，是將物質(zhì)波的概念和波動(dòng)方程相結(jié)合建立的二階偏微分方程，可描述微觀粒子的運(yùn)動(dòng)，每個(gè)微觀系統(tǒng)都有一個(gè)相應(yīng)的薛定諤方程式，通過解方程可得到波函數(shù)的具體形式以及對(duì)應(yīng)的能量，從而了解微觀系統(tǒng)的性質(zhì)。其僅適用于速度不太大的非相對(duì)論粒子，其中也沒有包含關(guān)于粒子自旋的描述。當(dāng)計(jì)及相對(duì)論效應(yīng)時(shí)，薛定諤方程由相對(duì)論量子力學(xué)方程所取代，其中自然包含了粒子的自旋。
薛定諤方程建立于 1926年。它是一個(gè)非相對(duì)論的波動(dòng)方程。它反映了描述微觀粒子的狀態(tài)隨時(shí)間變化的規(guī)律，它在量子力學(xué)中的地位相當(dāng)于牛頓定律對(duì)于經(jīng)典力學(xué)一樣，是量子力學(xué)的基本假設(shè)之一。設(shè)描述微觀粒子狀態(tài)的波函數(shù)為Ψ（r，t），質(zhì)量為m的微觀粒子在勢場V（r，t）中運(yùn)動(dòng)的薛定諤方程為
在給定初始條件和邊界條件以及波函數(shù)所滿足的單值、有限、連續(xù)的條件下，可解出波函數(shù)Ψ（r，t）。由此可計(jì)算粒子的分布概率和任何可能實(shí)驗(yàn)的平均值（期望值）。當(dāng)勢函數(shù)V不依賴于時(shí)間t時(shí)，粒子具有確定的能量，粒子的狀態(tài)稱為定態(tài)。定態(tài)時(shí)的波函數(shù)可寫成式中Ψ（r）稱為定態(tài)波函數(shù)，滿足定態(tài)薛定諤方程，這一方程在數(shù)學(xué)上稱為本征方程，式中E為本征值，是定態(tài)能量，Ψ（r）又稱為屬于本征值E的本征函數(shù)。
量子力學(xué)中求解粒子問題常歸結(jié)為解薛定諤方程或定態(tài)薛定諤方程。薛定諤方程揭示了微觀物理世界物質(zhì)運(yùn)動(dòng)的基本規(guī)律，被廣泛地用于原子物理、核物理和固體物理，對(duì)于原子、分子、核、固體等一系列問題中求解的結(jié)果都與實(shí)際符合得很好。
定態(tài)薛定諤方程直角坐標(biāo)系形式
定態(tài)薛定諤方程球坐標(biāo)系形式
1.2定態(tài)薛定諤方程
條件
V（r,t）=V(r)， 與t無關(guān)。
用分離變量法, 令Ψ=φ(r)f(t)，代入薛定諤方程，得兩個(gè)方程：
此稱定態(tài)薛定諤方程
整個(gè)定態(tài)波函數(shù)形式：
特點(diǎn)：
波函數(shù)由空間部分函數(shù)與時(shí)間部分函數(shù)相乘；
B．時(shí)間部分函數(shù)是確定的。
定態(tài)波函數(shù)幾率密度W與t無關(guān)，幾率分布不隨時(shí)間而變，因此稱為定態(tài)。
1.3本征方程、本征函數(shù)與本征值
算符： 本征方程：
λ：本征值，有多個(gè)，甚至無窮多個(gè)
ψλ：本征值為λ的本征函數(shù)，也有多個(gè)，甚至無窮多個(gè)，有時(shí)一個(gè)本征值對(duì)應(yīng)多個(gè)不同的本征函數(shù)，這稱為簡并。若一個(gè)本征值對(duì)應(yīng)的不同本征函數(shù)數(shù)目為N，則稱N重簡并。
1.4 定態(tài)情況下的薛定諤方程一般解
1、定態(tài)薛定諤方程或不含時(shí)的薛定諤方程是能量本征方程，E就稱為體系的能量本征值，而相應(yīng)的解稱為能量的本征函數(shù)。
2、當(dāng)不顯含時(shí)時(shí)，體系的能量是收恒量，可用分離變量。
3、解定態(tài)薛定諤方程，關(guān)鍵是寫出哈密頓量算符。
2. 利用矩陣法求解薛定諤方程
以在一維空間運(yùn)動(dòng)的粒子構(gòu)成的諧振子的體系為例。
該粒子的勢能是，是諧振子的角頻率，因此諧振子的哈密頓量為
。
當(dāng)時(shí)，諧振子的勢能變?yōu)闊o窮大，因此，粒子只能在有限的空間上運(yùn)動(dòng)，并且能量值譜是分立的。下面采用矩陣的方法，確定諧振子的能量分立值。
從運(yùn)動(dòng)方程出發(fā) （1）
而勢能 那么
又代入上式（1）得
即 （2）
在矩陣形式下，該方程可以寫為
含時(shí)坐標(biāo)矩陣元 （3）
對(duì)它求導(dǎo)，我們得到
代入上式后，有
（4）
其中 （5）
所以，除了當(dāng)或外，所有的坐標(biāo)矩陣元都等于零
當(dāng)時(shí)，由（5）式有
即 同理，
因此，只有變化時(shí)，才能得到頻率即 所以不為零的坐標(biāo)矩陣元為
根據(jù)定義[12-14]
對(duì)于存在的波函數(shù)，應(yīng)為實(shí)數(shù)，所有的矩陣元也為實(shí)數(shù)，由厄密算符的性質(zhì)得
為了計(jì)算坐標(biāo)的矩陣元，由對(duì)易關(guān)系
又 代入上式易得
寫為矩陣形式，有
根據(jù)矩陣的乘法規(guī)則，有
又，則有由前面的分析知，只有時(shí)，才存在矩陣元，代入上式，
從該方程我們可以得出
矩陣元不為零，但是當(dāng)時(shí)，矩陣元?jiǎng)t
即
又
依次類推，得出
最后，我們得到坐標(biāo)矩陣元不為零的表達(dá)式
又諧振子的能量可以用來表示，且，計(jì)算該能量得
其中，對(duì)于全部的1求和，只有當(dāng)參數(shù)時(shí)坐標(biāo)矩陣元不為零，因此得到
亦即
因此，諧振子的能級(jí)以為間隔，最低能級(jí)是
MATLAB仿真結(jié)果
線性諧振子的前六個(gè)本征函數(shù)
上圖為線性諧振子的前六個(gè)本征函數(shù)，圖中縱軸橫線表示具有相同能量的經(jīng)典線性諧振子的振動(dòng)范圍。
有限方勢阱前六個(gè)本征函數(shù)
上圖為有限方勢阱的前六個(gè)本征函數(shù)，圖中縱軸橫線表示具有相同能量的經(jīng)典線性諧振子的振動(dòng)范圍。
參考文獻(xiàn)：
1．周世勛，量子力學(xué)教程，北京-高等教育出版社，1979：38-42
2．曾謹(jǐn)言，量子力學(xué)，北京-科學(xué)出版社，1987：45-51
3. 周豐，定態(tài)薛定諤方程的計(jì)算機(jī)解法，武漢交通職業(yè)學(xué)院學(xué)報(bào)，2005，7（2），77-80
4. 封國林等，試用矩陣連分法數(shù)值求解薛定諤方程，江蘇農(nóng)學(xué)院學(xué)報(bào)，1996，17（4），103-108
5. 馬文淦 編著，計(jì)算物理學(xué)，科學(xué)出版社，2005，P196-201，244-250
6．王肇慶、佘守憲、蘇惠惠，諧振子薛定諤方程的簡單解法，大學(xué)物理，1996，15（8）：19-21
附錄：
程序運(yùn)行環(huán)境：MATLAB7.0
MATLAB源程序：
function f = schrodinger()
%% 對(duì)一些常數(shù)的定義
me = 9.10938188e-31;
eV = 1.60217646e-19;
h = 6.626068e-34;
hbar = 1.05457148e-34; % hbar=h/2/pi
% 定義寬度和格點(diǎn)數(shù)
a = 10e-9; % 長度
n = 128; % 分離數(shù)
z = linspace(-a/2,a/2,n); % 線性等分
dz = a/n; % 各點(diǎn)空間
%% 可能的矩陣
%為有限方勢阱====================================
%V0 = 0*eV;
%V = -V0*ones(n,1);
% 線性諧振子 ==================================
K = 1;
V0=1*eV;
V =V0+1/2*K*z'.^2;
pmatrix = spdiags(V,0,n,n); % 創(chuàng)建稀疏矩陣
%% 用薛定諤矩陣求波函數(shù)
vector = zeros(n,3);
vector(1:n,1) = -hbar^2/(2*me)/dz^2;
vector(1:n,2) = 2*hbar^2/(2*me)/dz^2;
vector(2:n,3) = -hbar^2/(2*me)/dz^2;
vmatrix = spdiags(vector,-1:1,n,n);
matrix = pmatrix+vmatrix;
eignum = 6; % 設(shè)置特征值個(gè)數(shù)
% 求薛定諤方程的特征值
[eigvector, eigvalue] = eigs(matrix,eignum,0); %求指定的幾個(gè)特征值
diag(eigvalue)/eV %矩陣對(duì)角元素提取、創(chuàng)建對(duì)角矩陣
for i = 1:eignum,
wavefunction = eigvector(:,i);
energy = eigvalue(i,i);
%將engivector常規(guī)化
wavefunction = wavefunction/sqrt(sum(abs(wavefunction.^2)*dz));
% 作圖
figure(1);
subplot(eignum/2,2,i),plot(z,wavefunction);%創(chuàng)建子圖、畫波函數(shù)圖
% figure(2);
% plot(z,energy);%畫能量線圖
end
程序運(yùn)行結(jié)果：
>> schrodinger()
Iteration 1: a few Ritz values of the 20-by-20 matrix:
0
0
0
0
0
0
Iteration 2: a few Ritz values of the 20-by-20 matrix:
1.0e+018 *
1.3147
1.5308
1.8346
2.2932
3.0648
4.6360
Iteration 3: a few Ritz values of the 20-by-20 matrix:
1.0e+018 *
1.3147
1.5308
1.8346
2.2932
3.0648
4.6360
ans =
4.7476
4.0774
3.4021
2.7218
2.0365
1.3463
【定態(tài)薛定諤方程的MATLAB求解(一)】相關(guān)文章：
MATLAB正交振幅調(diào)制解調(diào)仿真分析(一)12-26
基于MATLAB的正交振幅調(diào)制與解調(diào)仿真分析(一)03-07
LabVIEW與MATLAB混合編程03-19
基于matlab的節(jié)電模糊控制系統(tǒng)03-07
論文寫作格式要求解析03-14