摘要：此文介紹了一種使用MATLAB求解一維定態(tài)薛定諤方程的方法。利用充分格式進(jìn)行離散化，得出相應(yīng)的矩陣方程，用MATLAB求解本征值和本征函數(shù)。此方法簡(jiǎn)單可靠，可以處理各種時(shí)間無(wú)關(guān)的束縛態(tài)問(wèn)題。所用的程序附于文后。

---

宏觀物體遵循的是牛頓運(yùn)動(dòng)方程，給定初始條件以及受力條件，我們就可以求出任意時(shí)刻粒子的狀態(tài)。而在原子尺度上，所有粒子都表現(xiàn)出波的行為，粒子的狀態(tài)用波函數(shù)

來(lái)描述，描述微觀粒子運(yùn)動(dòng)的方程是由薛定諤于1925年提出薛定諤方程。微觀粒子的運(yùn)動(dòng)與其所處的勢(shì)場(chǎng)相關(guān)，當(dāng)勢(shì)場(chǎng)不隨時(shí)間變化時(shí)，稱(chēng)為定態(tài)，一維情況下的定態(tài)薛定諤方程為

其中，

表示粒子所處的勢(shì)場(chǎng)。由定態(tài)波函數(shù)可以得出總波函數(shù)

1926年， Born提出，波函數(shù)模的平方

代表在位置 ，時(shí)刻 尋找到電子的概率，這就是波函數(shù)的條統(tǒng)計(jì)解釋。可以由遵循微分方程的波函數(shù)表示，薛定諤波方程涉及空間坐標(biāo)和時(shí)間。對(duì)于一維情況，在時(shí)刻 時(shí)，在 和 之間的某處找到電子的概率由下式給出

時(shí)間無(wú)關(guān)薛定諤方程(1)是系統(tǒng)的能量本征方程。該特征值方程通常由一組特定的函數(shù)

和相應(yīng)的一組常數(shù) 滿足解的條件，被稱(chēng)為是哈密頓算符的本征函數(shù)和相應(yīng)的本征值。當(dāng)測(cè)量處于 狀態(tài)的系統(tǒng)的能量時(shí)，結(jié)果將始終為 。對(duì)于束縛態(tài)情況，必須有

一維定態(tài)薛定諤方程是二階微分方程，但是，能夠解析求解的情況屈指可數(shù)，如氫原子，諧振子，無(wú)限深勢(shì)阱等。隨著計(jì)算機(jī)技術(shù)的發(fā)展，我們可以數(shù)值上進(jìn)行求解。本文利用MATLAB軟件，使用矩陣方法求解束縛態(tài)的本征值問(wèn)題。

對(duì)于原子系統(tǒng)，以nm和eV的能量測(cè)量長(zhǎng)度更方便。我們可以使用縮放因子

因此我們可以將等式(1)寫(xiě)成

為了求解上述方程式，首先進(jìn)行離散化處理。將坐標(biāo)

離散化為 個(gè)點(diǎn)，用 來(lái)表示，若 ，則有 另外，還需要對(duì)方程式進(jìn)行差分格式處理，二階導(dǎo)數(shù)處理如下

因此，

的二階導(dǎo)數(shù)矩陣可以寫(xiě)成 矩陣

矩陣大小為 而不是 ，因?yàn)楹瘮?shù)的二階導(dǎo)數(shù)無(wú)法在終點(diǎn)進(jìn)行計(jì)算，即 和 。動(dòng)能矩陣為 ，勢(shì)能矩陣為對(duì)角矩陣，即 。則我們現(xiàn)在可以將哈密頓矩陣定義為 。用于生成哈密頓矩陣的代碼是% Make Second Derivative Matrix ---------------

因此，矩陣形式的薛定諤方程是

。MATLAB有內(nèi)置函數(shù)可以求解本征值問(wèn)題，其代碼為[e_funct, e_values] = eig(H)

其中e_funct是具有對(duì)應(yīng)于第n個(gè)本征函數(shù)的第n列的

矩陣，并且e_values是按遞增順序的N個(gè)本征值的列向量。一般可以求解出N-2個(gè)本征值和本征函數(shù)，但只有e_values的負(fù)值才有意義。為了獲得完整的特征向量，我們需要包括端點(diǎn)，其中 。

接下來(lái)舉一個(gè)例子，計(jì)算有限深方勢(shì)阱問(wèn)題。如圖所示是求解得到的有限深方勢(shì)阱的本征能量譜。

有限深方勢(shì)阱的本征能量譜

同時(shí)還可以求出本征函數(shù)以及幾率分布，如圖所示，

本征波函數(shù)以及幾率分布

我們還可以改變勢(shì)函數(shù)去計(jì)算各種各樣的定態(tài)束縛態(tài)問(wèn)題，也可以去計(jì)算已知解的問(wèn)題，以驗(yàn)證此方法的可靠性。

---

程序：

clear; clc; tic; num = 2001; % Number of data points (odd number) % Constants ----------------- hbar = 1.055e-34; e = 1.602e-19; m = 9.109e-31; eps0 = 8.854e-12; Ese = 1.6e-19; % Energy scaling factor Lse = 1e-9; % Length scaling factor Cse = -hbar^2/(2*m) / (Lse^2*Ese); % Schrodinger Eq constant % Potential well parameters U = zeros(num,1); U_matrix = zeros(num-2); % Potential Wells square well % Enter energies in eV and distances in nm xMin = -0.1; xMax = +0.1; x1 = 0.05; % 1/2 well width U1 = -400; % Depth of well (eV) x = linspace(xMin,xMax, num); for cn = 1 : numif abs(x(cn)) <= x1U(cn) = U1;end end s = sprintf('Potential Well: SQUARE'); % Graphics ----------------------- figure(1); set(gcf,'Name','Potential Energy','NumberTitle','off') plot(x,U,'LineWidth',3); axis([xMin-eps xMax min(U)-50 max(U)+50]);title(s); xlabel('x (nm)'); ylabel('energy (eV)'); grid on% Make potential energy matrix dx = (x(2)-x(1)); dx2 = dx^2; for cn =1:(num-2)U_matrix(cn,cn) = U(cn+1); end % Make Second Derivative Matrix off = ones(num-3,1); SD_matrix = (-2*eye(num-2) + diag(off,1) + diag(off,-1))/dx2; % Make KE Matrix K_matrix = Cse * SD_matrix; % Make Hamiltonian Matrix H_matrix = K_matrix + U_matrix; % Find Eignevalues E_n and Eigenfunctions psi_N [e_funct, e_values] = eig(H_matrix); % All Eigenvalues 1, 2 , ... n where E_N < 0 flag = 0; n = 1; while flag == 0E(n) = e_values(n,n);if E(n) > 0flag = 1;end % ifn = n + 1; end % while E(n-1) = []; n = n-2; % Corresponding Eigenfunctions 1, 2, ... ,n: Normalizing the wavefunction for cn = 1 : npsi(:,cn) = [0; e_funct(:,cn); 0]; area = simpson1d((psi(:,cn) .* psi(:,cn))',xMin,xMax);psi(:,cn) = psi(:,cn)/sqrt(area); % normalizeprob(:,cn) = psi(:,cn) .* psi(:,cn);if psi(5,cn) < 0psi(:,cn) = -psi(:,cn); end % curve starts positive end % for % Display eigenvalues in Command Window disp(' '); disp('========================= '); disp(' '); fprintf('No. bound states found = %0.0g n',n); disp(' '); disp('Quantum State / Eigenvalues En (eV)'); for cn = 1 : nfprintf(' %0.0f ',cn);fprintf(' %0.5g n',E(cn)); end disp(' ') disp(' ');% Plot energy spectrum xs(1) = xMin; xs(2) = xMax; figure(2); set(gcf,'Units','Normalized'); set(gcf,'Position',[0.5 0.1 0.4 0.6]); set(gcf,'Name','Energy Spectrum','NumberTitle','off') set(gcf,'color',[1 1 1]); set(gca,'fontSize',12); plot(x,U,'b','LineWidth',2); xlabel('position x (nm)','FontSize',12); ylabel('energy U, E_n (eV)','FontSize',12); h_title = title(s); set(h_title,'FontSize',12); hold on cnmax = length(E); for cn = 1 : cnmaxys(1) = E(cn);ys(2) = ys(1);plot(xs,ys,'r','LineWidth',2); end %for axis([xMin-eps xMax min(U)-50 max(U)+50]);% Plots first 5 wavefunctions & probability density functions if n < 6nMax = n; elsenMax = 5; end figure(3) clf set(gcf,'Units','Normalized'); set(gcf,'Position',[0.05 0.1 0.4 0.6]); set(gcf,'NumberTitle','off'); set(gcf,'Name','Eigenvectors & Prob. densities'); set(gcf,'Color',[1 1 1]); %nMax = 8; for cn = 1:nMaxsubplot(nMax,2,2*cn-1);y1 = psi(:,cn) ./ (max(psi(:,cn)-min(psi(:,cn))));y2 = 1 + 2 * U ./ (max(U) - min(U));plot(x,y1,'lineWidth',2)hold onplot(x,y2,'r','lineWidth',1)%plotyy(x,psi(:,cn),x,U);axis off%title('psi cn);title_m = ['psi n = ', num2str(cn)] ;title(title_m,'Fontsize',10);subplot(nMax,2,2*cn);y1 = prob(:,cn) ./ max(prob(:,cn));y2 = 1 + 2 * U ./ (max(U) - min(U));plot(x,y1,'lineWidth',2)hold onplot(x,y2,'r','lineWidth',1)title_m = ['psi^2 n = ', num2str(cn)] ;title(title_m,'Fontsize',10);axis off end toc

m文件simpson1d.m

function integral = simpson1d(f,a,b) % [1D] integration - Simpson's 1/3 rule % f function a = lower bound b = upper bound % Must have odd number of data points % Simpson's coefficients 1 4 2 4 ... 2 4 1 numS = length(f); % number of data points if mod(numS,2) == 1sc = 2*ones(numS,1);sc(2:2:numS-1) = 4;sc(1) = 1; sc(numS) = 1;h = (b-a)/(numS-1);integral = (h/3) * f * sc; else integral = 'Length of function must be an ODD number' end

參考資料：

http://www.physics.usyd.edu.au/teach_res/mp/mphome.htm

總結(jié)

以上是生活随笔為你收集整理的matlab向量的模_基于MATLAB使用矩阵方法求解一维定态薛定谔方程的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。