摘要：此文介紹了一種使用MATLAB求解一維定態薛定諤方程的方法。利用充分格式進行離散化，得出相應的矩陣方程，用MATLAB求解本征值和本征函數。此方法簡單可靠，可以處理各種時間無關的束縛態問題。所用的程序附于文后。

---

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

來描述，描述微觀粒子運動的方程是由薛定諤于1925年提出薛定諤方程。微觀粒子的運動與其所處的勢場相關，當勢場不隨時間變化時，稱為定態，一維情況下的定態薛定諤方程為

其中，

表示粒子所處的勢場。由定態波函數可以得出總波函數

1926年， Born提出，波函數模的平方

代表在位置 ，時刻 尋找到電子的概率，這就是波函數的條統計解釋?？梢杂勺裱⒎址匠痰牟ê瘮当硎?#xff0c;薛定諤波方程涉及空間坐標和時間。對于一維情況，在時刻 時，在 和 之間的某處找到電子的概率由下式給出

時間無關薛定諤方程(1)是系統的能量本征方程。該特征值方程通常由一組特定的函數

和相應的一組常數 滿足解的條件，被稱為是哈密頓算符的本征函數和相應的本征值。當測量處于 狀態的系統的能量時，結果將始終為 。對于束縛態情況，必須有

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

對于原子系統，以nm和eV的能量測量長度更方便。我們可以使用縮放因子

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

為了求解上述方程式，首先進行離散化處理。將坐標

離散化為 個點，用 來表示，若 ，則有 另外，還需要對方程式進行差分格式處理，二階導數處理如下

因此，

的二階導數矩陣可以寫成 矩陣

矩陣大小為 而不是 ，因為函數的二階導數無法在終點進行計算，即 和 。動能矩陣為 ，勢能矩陣為對角矩陣，即 。則我們現在可以將哈密頓矩陣定義為 。用于生成哈密頓矩陣的代碼是% Make Second Derivative Matrix ---------------

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

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

其中e_funct是具有對應于第n個本征函數的第n列的

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

接下來舉一個例子，計算有限深方勢阱問題。如圖所示是求解得到的有限深方勢阱的本征能量譜。

有限深方勢阱的本征能量譜

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

本征波函數以及幾率分布

我們還可以改變勢函數去計算各種各樣的定態束縛態問題，也可以去計算已知解的問題，以驗證此方法的可靠性。

---

程序：

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。