怎樣利用MATLAB求一維含時薛定諤方程的數值解？

01-07

寫一個和 @Narayan 不同的思路。既然是使用matlab，一定要用矩陣的思維！

首先，注意到這個問題的希爾伯特空間是無窮維的(空間上)，跟我先念:"無窮維問題先要化為有窮維"。無窮維問題怎麼化為有限維呢？差分。

首先我們確定問題的維數以及其他的參數。

L = 2*pi; % 範圍 N = 100; %有限維問題的維數 x = linspace(0,L,N)"; % 離散化坐標 dx = x(2) - x(1); %步長

我們在量子力學課上學到，算符一般都有其矩陣表示。

幾個比較原理性的算符，一階導數和二階導數算符的矩陣表示怎麼寫呢?

注意到 $Du approx frac{u(x+Delta x)-u(x-Delta x)}{2Delta x}$

我們如下構造導數算符 $D$ ：

D=(diag(ones((N-1),1),1)-diag(ones((N-1),1),-1))/(2*dx);

作為希爾伯特空間這樣一個特殊的空間，我們有一些約束條件的(空間自己的性質，所以應當表現在運算元自己的性質上)！因此我們有:

D(1,1) = 0; D(1,2) = 0; D(2,1) = 0; % So that f(0) = 0 D(N,N-1) = 0; D(N-1,N) = 0; D(N,N) = 0; % So that f(L) = 0

然後構造Laplacian運算元(exercise):

Lap = (-2*diag(ones(N,1),0) + diag(ones((N-1),1),1) ... + diag(ones((N-1),1),-1))/(dx^2); % Next modify Lap so that it is consistent with f(0) = f(L) = 0 Lap(1,1) = 0; Lap(1,2) = 0; Lap(2,1) = 0; % So that f(0) = 0 Lap(N,N-1) = 0; Lap(N-1,N) = 0; Lap(N,N) = 0;% So that f(L) = 0

現在作為演示，我們來試試這兩個運算元吧~

f = sin(x); % And try the following: Df = D*f; Lapf = Lap*f; plot(x,f,"b",x,Df,"r", x,Lapf,"g"); axis([0 5 -1.1 1.1]); % Optimized axis parameters % To display the matrix D on screen, simply type D and return ... D % Displays the matrix D in the workspace Lap % Displays the matrix Lap

Great!

下面先解定態玩玩先，注意，定態薛定諤方程是一個斯圖姆劉維爾本徵值問題:既然是一個本徵值問題，有了矩陣表示就直接用matlab超快的本徵系統求解器eig！

hbar = 1; m = 1; H = -(1/2)*(hbar^2/m)*Lap; % Solve for eigenvector matrix V and eigenvalue matrix E of H [V,E] = eig(H); % Plot lowest 3 eigenfunctions plot(x,V(:,3),"r",x,V(:,4),"b",x,V(:,5),"k"); shg; E % display eigenvalue matrix diag(E) % display a vector containing the eigenvalues

Narayan用了有限差分法繼續去做，我這裡則提供另一個方法，在勢不含時的情況下，演化是由演化算符 $U(t)= heta(t)e^{iHt}$ 控制的，注意在matlab里，一個矩陣的e指數使用expm命令構建。因此，如下代碼演示了一個波包傳播至50處的階梯然後散射和反射的過程。:

L = 100; % Interval Length N = 400; % No of points x = linspace(0,L,N)"; % Coordinate vector dx = x(2) - x(1); % Coordinate step % Parameters for making intial momentum space wavefunction phi(k) ko = 2; % Peak momentum a = 20; % Momentum width parameter dk = 2*pi/L; % Momentum step km=N*dk; % Momentum limit k=linspace(0,+km,N)"; % Momentum vector % Make psi(x,0) from Gaussian kspace wavefunction phi(k) using % fast fourier transform : phi = exp(-a*(k-ko).^2).*exp(-i*6*k.^2); % unnormalized phi(k) psi = ifft(phi); % multiplies phi by expikx and integrates vs. x psi = psi/sqrt(psi"*psi*dx); % normalize the psi(x,0) % Expectation value of energy; e.g. for the parameters % chosen above & = 2.062. avgE = phi"*0.5*diag(k.^2,0)*phi*dk/(phi"*phi*dk); % CHOOSE POTENTIAL U(X): Either U = 0 OR % U = step potential of height avgE that is located at x=L/2 %U = 0*heaviside(x-(L/2)); % free particle wave packet evolution U = avgE*heaviside(x-(L/2)); % scattering off step potential % Finite-difference representation of Laplacian and Hamiltonian e = ones(N,1); Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2; H = -(1/2)*Lap + spdiags(U,0,N,N); % Parameters for computing psi(x,T) at different times 0 &< T &< TF NT = 200; % No. of time steps TF = 29; T = linspace(0,TF,NT); % Time vector dT = T(2)-T(1); % Time step hbar = 1; % Time displacement operator E=exp(-iHdT/hbar) E = expm(-i*full(H)*dT/hbar); % time desplacement operator %*************************************************************** % Simulate rho(x,T) and plot for each T %*************************************************************** for t = 1:NT; % time index for loop % calculate probability density rho(x,T) psi = E*psi; % calculate new psi from old psi rho = conj(psi).*psi; % rho(x,T) plot(x,rho,"k"); % plot rho(x,T) vs. x axis([0 L 0 0.15]); % set x,y axis parameters for plotting xlabel("x [m]", "FontSize", 24); ylabel("probability density [1/m]","FontSize", 24); pause(0.05); % pause between each frame displayed end % Calculate Reflection probability R=0; for a=1:N/2; R=R+rho(a); end R=R*dx

參考文獻(代碼全部都是抄的):https://arxiv.org/pdf/0704.1622.pdf,MATLAB codes for teaching quantum physics.

（利用坐標表象下的場算符的矩陣表示也可以寫！為了讓沒學過場論初學者看懂，就不用場算符的矩陣表示這個詮釋了！）

簡單先寫一個Matlab解一維Schrodinger方程的例子吧：

首先 $ihbardfrac{partial}{partial t} psi=-dfrac{hbar^2}{2m}dfrac{partial^2}{partial x^2}psi+V(x)psi$ 從數學形式上看依然是熱傳導方程，

只不過是虛熱傳導方程。

用Matlab解偏微分方程的思路是，以差分代替微分。

設波函數為 $u(x,t)$ ,

$u_tapprox dfrac{u(x,t+Delta t)-u(x,t)}{Delta t}$ ;

$u_{xx}approx dfrac{u(x+Delta x,t)-2u(x,t)+u(x-Delta x,t)}{(Delta x)^2}$ ;

記 $x=iDelta x,t=jDelta t,dfrac{hbar^2}{2m}=a^2$ ,記虛數單位 $i o I$

計算轉化為對矩陣的計算，有遞推：

$dfrac{I}{Delta t}(u_{i,j+1}-u_{i,j})=-dfrac{a^2}{(Delta x)^2}(u_{i+1,j}-2u_{i,j}+u_{i-1,j})+V_{i,j}u_{i,j}$ ;

$u_{i,j+1}=a^2Idfrac{Delta t}{(Delta x)^2}(u_{i+1,j}+u_{i-1,j})-(V_{i,j}-2dfrac{a^2}{(Delta x)^2})I u_{i,j}$ ;

以下具體程序以勢壘為例（隧穿效應）：

x=1:220; v(x)=0; v(105:115)=0.18; %Potential k0=0.5; d=10; x0=40; %Parameter of Package B0=exp(k0*1i*x).*exp(-(x-x0).^2*log10(2)./d^2);%Gauss Package B(:,1)=B0."; A=diag(-2+2i+v)+diag(ones(219,1),1)+diag(ones(219,1),-1); for t=1:130 C(:,t+1)=4i*(AB(:,t)); B(:,t+1)=C(:,t+1)-B(:,t); plot(x,abs(B(:,t)).^2/norm(B(:,t)).^2,x,v/2); pause(0.1) end

部分截圖如下：