使用有限差分法求解无限长直的金属槽内的电位分布
实验内容：
实验源码：
实验截图

实验内容：
使用有限差分法求解无限长直的金属槽内的电位分布


实验源码：
%控制宽度
Length = 50;    %长度
Width = 50;     %宽度
Height = 1:1:50;

%控制轴的位置
Var_Xz = 25;
Var_Yz = 10;

Values1 = zeros(Width,Length);   %构造一个Width * Length的0矩阵
for m1 = 1 : Length              %将金属板盖设置为100V
    Values1(Width,m1)=100;
end

m=1;t=0;k=0;
Values2 = Values1;

while(m>0.001)
    m=0;
    k=k+1; %迭代次数
    
    %看做为网格，就是计算完成一个网格里的所有值与最初始的值对比误差
    for i=2:Width-1             %控制行；计算列值
        for j=2:Length-1        %控制列；计算行值    
            Values2(i,j)= (Values1(i,j+1) + Values1(i+1,j) + Values1(i-1,j) + Values1(i,j-1)) / 4;
            t=(Values2(i,j) - Values1(i,j));    % v2 与 v1 偏差
            if(t>m)                             %精度比较控制
                m=t;
            end 
        end
    end
    Values1=Values2;
end

%xoy面等势线截面
subplot(2,3,[1,2]);
contour(Values2,40);            %矩阵的等高线图
xlabel("x轴");
ylabel("y轴");
title(['1.xoy面等势线截面图','迭代次数：',num2str(k),'次'],'fontsize',10);

%xoy面电势分布
subplot(2,3,3);
face_xoy = surface(Values1);
face_xoy.FaceColor = 'flat';
face_xoy.EdgeColor = 'none';
axis([0,50,0,50]);
xlabel("x轴");
ylabel("y轴");
title("2.xoy面电势分布图",'fontsize',10);
%---------------------------------------------

%yoz面电势分布--当 x = Var_Xz
subplot(2,3,4);
yozValues = repmat(Values1(:,Var_Xz),1,50);
face_yoz = surface(yozValues);
face_yoz.EdgeColor = 'none';
axis([0,50,0,50]);
xlabel("z轴");
ylabel("y轴");
title(["3.yoz面电势分布-x = ",num2str(Var_Xz)],'fontsize',10);

%xoz面电势分布--当 y = Var_Yz
subplot(2,3,5);
xozValues = repmat(Values1(Var_Yz,:).',1,50);
face_xoz = surface(xozValues);
face_xoz.EdgeColor = 'none';
axis([0,50,0,50]);
xlabel("z轴");
ylabel("x轴");
title(["4.xoz面电势分布-y = ",num2str(Var_Yz)],'fontsize',10);

%截取一部分的立方体电势图
subplot(2,3,6);
[X, Y] = meshgrid(0:1:50, 0:1:50);
Z = zeros(size(X));
CO = Z; % blue
CO1 = Z * 50;

%xoy两面
xy1 = surf(X,Y,Z,Values1);hold on;
xy2 = surf(X,Y,Z+50,Values1);hold on;
%xoz两面
xz1 = surf(Y,Z,X,CO);hold on;
xz2 = surf(Y,Z+50,X,CO-10);hold on;
%yoz面
yz1 = surf(Z,Y,X,CO);hold on;
yz2 = surf(Z+50,Y,X,CO);hold on;

xz1.EdgeColor = 'none';
xz2.EdgeColor = 'none';
yz1.EdgeColor = 'none';
yz2.EdgeColor = 'none';
xy1.EdgeColor = 'none';
xy2.EdgeColor = 'none';

xlabel("x轴");
ylabel("y轴");
zlabel("z轴");
axis([0,60,0,60,0,60]);
title("5.带盖金属3D电势分布图");

实验截图


有需要的可以一起交流学习、觉得有帮助的话，点个赞再走啦、谢谢啦

在前面的笔记里孤光一点萤：数值常微分方程-欧拉法与龙格-库塔法​zhuanlan.zhihu.com
整理了常微分方程的一些数值解法。类似的方法可以拓展到解偏微分方程的问题，这里整理有限差分法的相关笔记。
计算机无法处理连续介质的无限多点，因此只能通过网格划分的方式，将带求解区域内连续分布的函数进行离散化，这样只用求解有限多的格点即可。同样地，在计算机里无法取极限求导数，只能通过离散化之后，用有限差商代替导数。由此，可以将带求解的偏微分方程转化为一组离散格点上的差分方程。解这组方程的问题，最终又归化为一个矩阵的求解问题，而矩阵求解这种“体力活”是计算机所擅长的。
由于我们使用划分的离散网格来近似连续的待求解空间，因此网格划分越密，求解越精确。同时，网格划分越密集，意味着最后得到的矩阵方程维度越高，计算量会非常庞大。这时候常用的高斯消元法(一个高斯消元的例子)往往会吃不消，需要使用松弛迭代法等。
首先，还是从单变量的方程开始。考虑函数
 ，其中 
 。将区间 
 分为N份，每一份就是步长 
 ，则离散化后的格点坐标为： 
 。对于任意格点 
 ：
一阶向前差分
一阶向后差分
一阶中心差分
由此可得：
一阶向前差商
一阶向后差商
一阶中心差商
二阶中心差商
现在拓展到二维情况：
例如对于方程
 ，待求解区域为正方形域 
将待求解区域划分为矩形网格，每一个小矩形的长和宽(步长)分别为
 和 
 ，使用中心差商，可得：
对于这里的正方形区域，自然而然可以选取
 ，代入方程，得：
以此类推，可以得到对于点
 ，有：
这便是点
 满足的差分方程，由于该方程包含点 
 及其上下左右一共5个点，也称为“五点差分格式”。
五点差分格式的截断误差为
 阶。如果把左上、左下、右上、右下四个点也用上，为“九点差分格式”，误差可以降低到 
 阶。
为简便起见，这里以拉普拉斯方程
 为例，则五点差分格式简化为：
在边界
 上， 
在划分网格时候要注意，如果划分网格为
 个正方形格子，那么待求解的点有 
 个，因为区域边界上的点(即边界条件)是给定的，不用求。因此，给网格编号时候，从0开始编号再合适不过了，如下：
考虑第一个点
 ，根据五点差分格式，它的值由上下左右四个点确定，很显然，上、左分别为边界点 
 和 
 ，右、下分为为点 
 和 
 ，写为差分格式为：
第一列其他点同理，以此类推，如果将网格第
 列上每个点的解记为：
 ，则最终可以写为：
其中
 为 
 阶单位矩阵，
这里出现了一个迭代方程，想解
 需要 
 ，解 
 需要 
 ......
然后把这一系列矩阵方程写成一个更大的矩阵方程：
把整个网格的每一列纵向连接为一个列向量：
则有：
 ，其中：
 ， 
终于，我们把要求解的问题化简为了矩阵方程
 的形式，可以方便地使用计算机求解，这里A是系数矩阵，列向量 
 便是问题的解。
对于上面的拉普拉斯方程，用Python求解：
为了方便，求解时直接使用了线性代数库
 ，这样计算的是矩阵方程的精确解，耗时有点长。可以用松弛迭代法来求数值解，有空再更。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
N=79 # 矩阵维数
Imat = -np.eye(N)/4
def Gmatrix(N): # 生成G矩阵
mtr = np.eye(N)
for i in range(0, N - 1):
mtr[i][i + 1] = -1 / 4
for i in range(1, N):
mtr[i][i - 1] = -1 / 4
return mtr
def zero_mat(N,m): # 生成m个横向连接N阶零矩阵,m>=1
zero_matrix = np.zeros((N,N))
for i in range(m-1):
zero_matrix = np.block([zero_matrix,np.zeros((N,N))])
return zero_matrix
def zero_arr(N, m):  # 生成m个纵向连接N阶列向量,m>=1
zero_array = np.zeros((N,1))
for i in range(m-1):
zero_array = np.block([[zero_array],[np.zeros((N,1))]])
return zero_array
def Kmatrix(N,Gmat,Imat):
Kmar_list = []
for i in range(N):
if i == 0:
Kmat_line = np.block([np.block([Gmat, Imat]),zero_mat(N,N-2)])
Kmar_list.append(Kmat_line)
elif i == 1:
Kmat_line = np.block([Imat, Kmar_list[-1][:,:-N]])
Kmar_list.append(Kmat_line)
else:
Kmat_line = np.block([np.zeros((N,N)), Kmar_list[-1][:,:-N]])
Kmar_list.append(Kmat_line)
Kmat = Kmar_list[0]
for j in range(N-1):
Kmat = np.block([[Kmat],[Kmar_list[j+1]]])
return Kmat
def Barray(N): # 方程右边
barr = np.ones((N,1))/4
barr = np.block([[barr],[zero_arr(N,N-2)]])
barr = np.block([[barr],[np.ones((N,1))/4]])
return barr
Gmat = Gmatrix(N)
Kmat = Kmatrix(N,Gmat,Imat)
Barr = Barray(N)
Phi = np.linalg.solve(Kmat,Barr)
def trans(N,Phi):
# 把线性方程组解得的列矢量分割，从左到右排列转化为方矩阵，方便可视化
phi = Phi[0:N]
for i in range(N-1):
j = i+1
phi = np.block([phi,Phi[(j*N):((j+1)*N)]])
return phi
Phi = trans(N,Phi)
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
X,Y = np.meshgrid(x,y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X,Y,Phi,rstride=1,cstride=1,cmap=plt.get_cmap('rainbow'))
ax.contourf(X,Y,Phi,zdir='z',offset=1,cmap='rainbow')
plt.title('Difference methods for Laplace equation',fontsize='12')
plt.show()
作图如下：
通过控制作图参数里的offset，可以调整等高线的位置：
差分法求解的误差主要来自于差分公式对偏微分的近似替代，为了提高求解精度，可以把网格划分得更密，同时计算量也会显著增大。
Reference：
1、马文淦，《计算物理学》
2、李荣华，刘播，《微分方程数值解法》

	

一维势阱是量子力学中最简单最基本的模型，在量子力学教学和科研过程中起着基础性作用[1]。其理论结果在许多实际系统中也得到了很好的应用，比如在低维量子系统(如量子点、量子面、量子线等)中[2]。但作为量子力学中的基本方程，薛定谔方程[3]的求解却并不简单，在《量子力学》课程的教学过程中，除了一维无限深势阱、一维谐振子势等特殊一维情形有解析解外，一般情形很难求解。即使对于处在一维有限深势阱中的运动的粒子，当其处于束缚态时，由于确定其能级的是超越方法，无法具体给出它们的能级的解析表达式和归一化波函数[4]。
薛定谔方程的求解可分为解析法和数值法两种，解析法主要有WKB法、变分法等，数值法有打靶法[5]和有限元(FEM)法等。作为一种有益的探索，本文以处理多层光波导问题的转移矩阵(TM)法为算法，求解粒子在不同一维有限深势阱中运动的能量本征值，并利用Matlab作为计算工具，绘制出对应的波函数图像。采用TM方法计算，通过不同势阱的设置，可以得到定性的准解析规律性结论，为解释实验现象提供重要的理论依据。而且，在《量子力学》课程的学习中，学生可借助不同势阱的计算分析将理论学习进行拓展，更直观更透彻地理解一维势阱相关方面的知识点，加深对相关物理概念的理解。
1 转移矩阵法
1.1 算法应用依据
转移矩阵法起源于光学，用于计算多层薄膜介质的反射率和透射率。在光波导技术中，转移矩阵法用来建立不同介质层间的场分布联系。而经过对比发现：不含时的薛定谔方程与平板波导的波动方程在形式上是一致的；同时一维势阱模型与平板波导模型相类似，因而转移矩阵法也能够用于解决一维有限深势阱问题[6]。平板波导的波动方程为[7](1)其中k0=2π/λ，是光在真空中的传播常数，Ey为y方向的电场，β为电磁场沿z方向的传播常数，nj对应于波导中不同层的折射率，j=0,1,2。不含时的一维薛定谔方程形式为其中V(x)表示势能，m为粒子质量，E和ψ分别代表粒子的能量本征值和本征函数，上式经过变形，有：(2)对比方程(2)与方程(1)即可发现：这两个方程在形式上是一致的，方程(2)中的本征函数ψ即对应方程(1)中的电场强度Ey，因此两个方程可采用相同的方法求解。按照波动光学理论，光波是电磁横波，光波在空间任意位置的电磁场强度及所在介质性能之间的联系都是通过Maxwell方程组和物质方程建立，电磁波在两种介质形成的界面上反射和透射时的振幅反射因数和透射因数均可以由菲涅耳公式确定。若是一个多界面的薄膜系统，则根据以下两点[8]：(1)在每一界面处运用电磁场边界条件，将同一界面两侧的场分布联系起来；(2)利用与电磁场传播相伴随的相位差，将同一膜层上下两界面内侧的场分布联系起来，可将多界面系统看作是入射介质与薄膜、基底形成的等效介质之间的界面。将以上思路应用到任何一个复杂的薄膜系统，此时系统中光反射率和透射率问题，都可以通过其等效界面对应的等效介质进行计算。等效介质的光谱特性可以采用一个特征矩阵来表示，该矩阵也称之为转移矩阵。非对称平板波导的模型如图1所示。通常来说折射率有n1>n0>n2，导波光可被约束在导波层中进行传播。将导波光类比为势场中运动的粒子，粒子会被约束在势能较低处运动，则导波层可看作势阱中势能低的地方，衬底层和覆盖层可看作势能相对较高的地方。从这一角度看，势阱中粒子运动的模型和波导中导波光的传播模型是一致的。基于以上分析，可将转移矩阵法应用到不含时薛定谔方程的求解。而在实际应用中，例如超晶格量子阱，虽然由能隙不同的材料组成可以形成复杂的势阱，但系统的基本特征完全可以通过有限深势阱问题获得很好的理解[2]。图1 非对称平板波导模型
1.2 建模与求解
对方程(1)的求解，首先选取两个特解E1(x)和E2(x)，使之满足E1(0)=E′2(0)=1和E′1(0)=E2(0)=0。则方程(1)的通解可设为E1(x)和E2(x)的线性叠加，即Ey(x)=C1E1(x)+C2E2(x)(3)在区间(0，d)界面上，根据(3)式可给出场分布及其导数，建立0-d界面上的转移关系：(4)M称为区间(0，d)的转移矩阵。定义若选取两个特解分别为cos(κx)和则M为：(5)可以看出，转移矩阵仅和区间内的折射率分布、传播常数以及d有关。借助转移矩阵M，界面x = 0和x = d处的电磁场分布建立起联系。而利用矩阵逆阵的概念，可得反向传递关系，有：(6)由于矩阵M及其逆阵是互为逆矩阵的关系，通常两者都被称为转移矩阵。转移矩阵法理论上可以求解任意的一维势阱，一维势阱的基本型就是一维方势阱或一维类方势阱。对于基本的一维方势阱，势能分布如图2所示。图2 一维方势阱对应束缚态的波函数为[7](7)其中(8)由转移矩阵理论，波函数ψ(x)及其一阶导数ψ′(x)应满足矩阵方程(9)根据束缚态波函数以及一维方势阱势能分布，可以得到波函数边界条件(10)将边界条件代入矩阵方程并进行一些简单的变换，就可得到方程(11)求解方程(11)即可获得粒子的能量E和相应的波函数ψ(x)。图3 一维任意势阱根据以上算法，可将一维方势阱推广到任意的一维任意势阱。对于一维任意势阱，其势能分布如图3所示。在x<0和x>d的区域，波函数已衰减到足够小，势场的变化带来的影响微乎其微，因此可在这两点截断，考虑x<0时，V(x)=V0，x>d时，V(x)=V2(V0，V2均为有限值)。在0<x<d的范围内，波函数为震荡形式。由于转移矩阵的建立需要依靠边界条件，而势阱内部的势能随x的变化而变化，因此需用到数学中的微元处理。如图3所示，在0<x<d的范围内将势能曲线化为n个相等宽度h的小区域，在每一个小区域内，势能视作一定值，由具体的势能函数确定。则每一个小区域都对应有一转移矩阵Mi，i=0，1，2，3，…，n。i=0,1,2，…，n(12)于是有(13)若在0<x<d中存在部分区域V>E，即势场能量高于粒子能量，此时可定义κi=iαi(14)利用恒等式sin(ix)=isinh(x),cos(ix)=cosh(x)则无需改变运算规则，同样可利用转移矩阵方法运算。可以发现计算过程中微元细分得越多，得到的能量本征值和波函数就越精确。
2 四种一维有限深势阱计算结果
图4 一维有限深势阱中波函数、位置分布概率图像绘制流程图图4给出了采用Matlab计算的流程图。对于不同的有限深势阱，在势阱形状参数设定后，可根据方程(11)求解粒子能量E值，如果E值满足束缚态条件，则继续运行得到势阱中的波函数和粒子的位置分布概率图示。
2.1 一维对称有限深势阱
图5(a)为一维对称有限深方势阱，图中势阱为2nm宽，两边势垒高度为2eV，势阱内以不同颜色的线从低到高分别来表示n=1,2,3所对应的基态能量和激发态能量。图5(b)和图5(c)则分别给出了相对应能级的波函数和位置分布概率。在实际物理问题中，微观粒子可以被势阱所束缚，但势阱并不是无限深，而是几个电子伏的有限深势阱。从图5中结果可以看出：有限深势阱中粒子的能量仍然是离散化的；在基态情况下，微观粒子最可能的位置仍然是位于势阱正中。此结果表明，无限深势阱虽是一个理想模型，但在一定条件下，很多系统都可以抽象为无限深势阱问题来处理。但有限深势阱结果和无限深势阱结果仍略有不同。由于微观粒子并不能穿透无限深势阱，因此在无限深势阱中粒子波函数表现为在边界上截断的正弦波，而有限深势阱中，无论是处于基态还是激发态，粒子均有一定几率穿出势阱，粒子波函数延续超出势阱边界，而且随着粒子能量的增加，穿出势阱的概率逐渐增大。
2.2 一维不对称有限深势阱
在三层材料体系中，若是左右两边生长的薄层材料不一样，则为不对称有限深势阱，如常见的非对称平板波导结构。图6(a)给出了不对称有限深势阱的典型结构，为便于与图5结果对比，势阱仍设置为2nm宽，但左势垒高度降为1eV，而右势垒高度保持为2eV。图6(b)给出n=1,3,5时的粒子波函数结果，可以看出：当粒子能量较低时，如n=1,3时，此时不对称势阱对粒子的波函数并无太大影响；而当粒子能量较高，如n=5时，粒子能量高于左势垒，此时的波函数连续进入左势垒，类似波导结构中的“辐射模”。而从图6(c)来看，不对称势阱对粒子的位置分布概率也有影响，n=1基态时位置分布概率不对称情形并不明显，但随着能量的升高，如n=3时，位置分布概率的不对称性明显增加，粒子进入低势垒侧的几率增大，而对于n=5，粒子有更高的概率存在于低势垒中。
2.3 一维类三角势阱
超晶格量子阱物理自20世纪70年代以来得到了长足的发展，不同形式的载流子的运动规律及在外场作用下的输运问题得到较好的研究[9]。在半导体器件中，单个异质界面，如MOS结构中Si/SiO2，GaAs/AlGaAs界面附近，载流子被限制在一很窄的势阱中，通常将这类势阱称为类三角量子阱，这是在电子器件问题中最常遇到的情形[10]，同样可以采用TM理论来讨论类三角势阱中的波函数分布。图7 一维类三角势阱结果演示
(a) 势阱；(b) 波函数；(c) 位置分布概率图7给出了一维类三角势阱的形式、波函数与位置分布概率的演示，给出的能级分别为n=1,2,3。从图7(a)来看，此时势阱中能级分布并不是均匀的，随能级升高，能级间的能量间隔逐渐减小。在图7(b)中，粒子的波函数在势阱左端全部为零，这是因为在实际的Si/SiO2界面，势垒高度可到3eV，波函数向势垒区渗入的影响完全可以忽略，效果类似无限深势阱中的“刚性壁”。而从图7(c)来看，随能级升高，粒子的分布概率呈准周期性的振荡，振荡周期越来越大，振荡幅度也越来越大，粒子很容易穿出右侧势阱而到达势阱外。对比分析类三角势阱和方势阱，可以发现：方势阱的底部水平，因此能量一定时德布罗意波长是个常数，对应波函数相邻两个节点的间距相等，这满足正弦函数的数学表征[9]；而对于类三角势阱，势阱底部从左到右上升，粒子的能量与势能之差减小，因此能量一定时德布罗意波长从左到右不断增大，波函数相邻两个节点的间距越来越大，振荡周期就越来越大。
2.4 一维有限对称双阶梯势阱
在半导体多层材料体系中，不同材料间的缓冲层设计可用来减少材料间的晶格失配，降低材料的表面、界面复合，形成的双面异质结的能带结构，则对应图8(a)中的对称双阶梯势阱，也可以将它视为一个简化的双势阱模型[11]。图8 一维有限双阶梯势阱结果演示
(a) 势阱；(b) 波函数；(c) 位置分布概率图8(b)中给出了n=1,3,5时的波函数演示结果。从图(b)中可以看出，在对称势阱中，波函数的分布始终是对称的。当n=5时，粒子的相应能量高于较低势阱，此时波函数连续进入了中间势阱。而从图8(c)的位置分布概率来看，此时粒子有更高的概率进入在中间势阱中。
3 总结
本文将处理光导波问题的TM法成功运用到求解一维有限深势阱的不含时薛定谔方程中，得到了粒子在不同势阱内运动的能量本征值，并给出了粒子波函数和位置分布概率的图像，直观地反映出粒子在势阱内的分布情况。借助Matlab强大的计算功能，可根据实际的物理问题和物理图像，对四种不同的一维有限深势阱进行设置，直观地显示不同势阱中的计算结果，这对学生学习《量子力学》是一个很好的辅助应用，能够帮助学生更深入地理解粒子在势场内的位置分布概率、运动行为等，加深对量子知识的理解。参考文献
[1] 曾谨言. 量子力学(卷1)[M]. 北京: 科学出版社, 2013.
[2] 夏吾吉, 张林. 量子力学中无限深势阱问题的教学研究[J]. 大学物理, 2015, 34(2): 35-40.
XIA W J, ZHANG L. Research on infinite potential well in quantum mechanics teaching[J]. College Physics, 2015, 34(2): 35-40. (in Chinese)
[3] 李爱君, 周战荣, 朱海飞. 大学物理课程中“波函数 薛定谔方程”的教学探讨[J]. 物理与工程, 2018, 28(4): 78-82, 87.
LI A J, ZHOU Z R, ZHU H F. Discussions on the teaching of “Wave function and Schrodinger equation” in college physics[J]. Physics and Engineering, 2018, 28(4): 78-82, 87. (in Chinese)
[4] 尹建武, 张勇和, 李大农. 方势阱中束缚态粒子能级的数值方法和波函数的图示[J]. 大学物理, 2008, 27(3): 33-36.
YIN J W, ZHANG Y H, LI D N. The numerical methods and graphical displays of the particles bound in one-dimensional square potential wells[J]. College Physics, 2008, 27(3): 33-36. (in Chinese)
[5] 徐中辉, 钟阳万, 肖庆生, 等. 用打靶法求解一维薛定谔方程的定态解[J]. 西南民族大学学报(自然科学版), 2009, 35(6): 1143-1146.
XU Z H, ZHONG Y W, XIAO Q S, et al. The steady-state solution of one-dimensional Schrodinger equation with the shooting method[J]. Journal of Southwest University for Naitonalities (Natural Science Edition), 2009, 35(6): 1143-1146. (in Chinese)
[6] 王忆锋, 唐利斌.利用转移矩阵和MATLAB求解一维薛定谔方程的一种简捷方法[J]. 红外技术, 2010, 32(3): 177-180.
WANG Y F, TANG L B. A simple method for solution of one-dimensional Schrodinger equation with transfer matrix and Matlab[J]. Infrared Technology, 2010, 32(3): 177-180. (in Chinese)
[7] 曹庄琪.导波光学[M]. 北京: 科学出版社, 2007.
[8] 卢进军, 刘卫国. 光学薄膜技术[M]. 西安: 西北工业大学出版社, 2005.
[9] 罗强, 姜玉梅, 韩玖荣,等.一维无限深梯形势阱中微观粒子波函数和能级的半解析解[J]. 大学物理，2014, 33(6): 55-60.
LUO Q, JIANG Y M, HAN J R, et al. The semi-analytical solution of wave functions and energy levels of a micro-particle in one-dimensional infinite trapezoid potential well[J]. College Physics, 2014, 33(6): 55-60. (in Chinese)
[10] 叶良修. 小尺寸半导体器件的蒙特卡罗模拟[M]. 北京: 科学出版社, 1997.
[11] 梁霞霞, 李冀光, 董晨钟.一维双势阱问题的理论研究[J].原子与分子物理学报, 2008, 25(5): 1039-1041.
LIANG X X, LI J G, DONG C Z. Investigation of one-dimension double-potential well[J]. Journal of Atomic and Molecular Physics, 2008, 25(5): 1039-1041. (in Chinese)
基金项目: 武汉理工大学教学研究项目(w2016002和w2015053)资助。作者简介: 陈凤翔，女，武汉理工大学教授，主要从事物理教学科研工作，研究方程为光伏科学与工程，phonixchen79@whut.edu.cn。引文格式: 陈凤翔，曹功辉，汪礼胜. 一维有限深势阱的转移矩阵法求解[J]. 物理与工程，2020，30(2)：37-43,48.
END
更多精彩文章请点击下面“蓝字”标题查看：
读后感：教育中的现实和远方
王青教授：昨晚(6月9日)，清华电动力学期末考试
朱邦芬院士：“减负”误区及我国科学教育面临的挑战
关于新冠肺炎疫情防控期间物理类课程线上教学的调查报告
乐永康：新冠肺炎疫情防控下美国物理实验教学及中美情况对比
教育部高等学校大学物理课程教学指导委员会关于推进在线物理教育教学研究的工作
王青，郭应寿：清华大学《费曼物理学II》和《电动力学》混合式线上教学实践
《物理与工程》2020年第1期目录
《物理与工程》2019年第6期目录
《物理与工程》2019年第5期目录
《物理与工程》2019年第4期目录
《物理与工程》2019年第3期目录
顾牡：对于重新制定的《非物理类理工学科大学物理课程教学基本要求》的认识和体会
黄昆：永远怀念吴大猷老师
黄昆：我的研究生涯
朱邦芬院士：从基础科学班到清华学堂物理班
朱邦芬院士：对培养一流拔尖创新人才的思考
李学潜教授：物理是一种文化
李学潜教授：如何帮助物理系学生迈过从高三到大一这个坎
穆良柱：什么是物理及物理文化？
穆良柱：什么是ETA物理认知模型
穆良柱：什么是ETA物理教学法
吴国祯教授：我的国外研究生经历印象——应清华大学物理系“基科班20年·学堂班10年纪念活动”而写
陈佳洱，赵凯华，王殖东：面向21世纪，急待重建我国的工科物理教育
王亚愚教授：清华物理系本科人才培养理念与实践
葛惟昆教授：关于中外人才培养的几点思考
王青教授：小班教学与翻转课堂：《费曼物理学Ⅱ》的10年教学实践——纪念费曼先生百年诞辰
安宇教授：为什么传统的课堂讲授模式需要改变
安宇教授：其实教学就是积累的过程
刘玉鑫教授：关于本科生物理基础课程教学和教材编著的一些思考
沈乾若：重创理科教育的美加课程改革
Henderson C：美国研究基金支持下的物理教育研究及其对高等物理教育的影响
《物理与工程》期刊是专注于物理教育教学研究的学术期刊，是中国科技核心期刊，1981年创刊，欢迎踊跃投稿，期刊投审稿采编平台：
http://gkwl.cbpt.cnki.net
欢迎关注
《物理与工程》微信公众号

在前面的笔记里孤光一点萤：数值常微分方程-欧拉法与龙格-库塔法​zhuanlan.zhihu.com
整理了常微分方程的一些数值解法。类似的方法可以拓展到解偏微分方程的问题，这里整理有限差分法的相关笔记。
计算机无法处理连续介质的无限多点，因此只能通过网格划分的方式，将带求解区域内连续分布的函数进行离散化，这样只用求解有限多的格点即可。同样地，在计算机里无法取极限求导数，只能通过离散化之后，用有限差商代替导数。由此，可以将带求解的偏微分方程转化为一组离散格点上的差分方程。解这组方程的问题，最终又归化为一个矩阵的求解问题，而矩阵求解这种“体力活”是计算机所擅长的。
由于我们使用划分的离散网格来近似连续的待求解空间，因此网格划分越密，求解越精确。同时，网格划分越密集，意味着最后得到的矩阵方程维度越高，计算量会非常庞大。这时候常用的高斯消元法(一个高斯消元的例子)往往会吃不消，需要使用松弛迭代法等。
首先，还是从单变量的方程开始。考虑函数
 ，其中 
 。将区间 
 分为N份，每一份就是步长 
 ，则离散化后的格点坐标为： 
 。对于任意格点 
 ：
一阶向前差分
一阶向后差分
一阶中心差分
由此可得：
一阶向前差商
一阶向后差商
一阶中心差商
二阶中心差商
现在拓展到二维情况：
例如对于方程
 ，待求解区域为正方形域 
将待求解区域划分为矩形网格，每一个小矩形的长和宽(步长)分别为
 和 
 ，使用中心差商，可得：
对于这里的正方形区域，自然而然可以选取
 ，代入方程，得：
以此类推，可以得到对于点
 ，有：
这便是点
 满足的差分方程，由于该方程包含点 
 及其上下左右一共5个点，也称为“五点差分格式”。
五点差分格式的截断误差为
 阶。如果把左上、左下、右上、右下四个点也用上，为“九点差分格式”，误差可以降低到 
 阶。
为简便起见，这里以拉普拉斯方程
 为例，则五点差分格式简化为：
在边界
 上， 
在划分网格时候要注意，如果划分网格为
 个正方形格子，那么待求解的点有 
 个，因为区域边界上的点(即边界条件)是给定的，不用求。因此，给网格编号时候，从0开始编号再合适不过了，如下：
考虑第一个点
 ，根据五点差分格式，它的值由上下左右四个点确定，很显然，上、左分别为边界点 
 和 
 ，右、下分为为点 
 和 
 ，写为差分格式为：
第一列其他点同理，以此类推，如果将网格第
 列上每个点的解记为：
 ，则最终可以写为：
其中
 为 
 阶单位矩阵，
这里出现了一个迭代方程，想解
 需要 
 ，解 
 需要 
 ......
然后把这一系列矩阵方程写成一个更大的矩阵方程：
把整个网格的每一列纵向连接为一个列向量：
则有：
 ，其中：
 ， 
终于，我们把要求解的问题化简为了矩阵方程
 的形式，可以方便地使用计算机求解，这里A是系数矩阵，列向量 
 便是问题的解。
对于上面的拉普拉斯方程，用Python求解：
为了方便，求解时直接使用了线性代数库
 ，这样计算的是矩阵方程的精确解，耗时有点长。可以用松弛迭代法来求数值解，有空再更。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
N=79 # 矩阵维数
Imat = -np.eye(N)/4
def Gmatrix(N): # 生成G矩阵
mtr = np.eye(N)
for i in range(0, N - 1):
mtr[i][i + 1] = -1 / 4
for i in range(1, N):
mtr[i][i - 1] = -1 / 4
return mtr
def zero_mat(N,m): # 生成m个横向连接N阶零矩阵,m>=1
zero_matrix = np.zeros((N,N))
for i in range(m-1):
zero_matrix = np.block([zero_matrix,np.zeros((N,N))])
return zero_matrix
def zero_arr(N, m):  # 生成m个纵向连接N阶列向量,m>=1
zero_array = np.zeros((N,1))
for i in range(m-1):
zero_array = np.block([[zero_array],[np.zeros((N,1))]])
return zero_array
def Kmatrix(N,Gmat,Imat):
Kmar_list = []
for i in range(N):
if i == 0:
Kmat_line = np.block([np.block([Gmat, Imat]),zero_mat(N,N-2)])
Kmar_list.append(Kmat_line)
elif i == 1:
Kmat_line = np.block([Imat, Kmar_list[-1][:,:-N]])
Kmar_list.append(Kmat_line)
else:
Kmat_line = np.block([np.zeros((N,N)), Kmar_list[-1][:,:-N]])
Kmar_list.append(Kmat_line)
Kmat = Kmar_list[0]
for j in range(N-1):
Kmat = np.block([[Kmat],[Kmar_list[j+1]]])
return Kmat
def Barray(N): # 方程右边
barr = np.ones((N,1))/4
barr = np.block([[barr],[zero_arr(N,N-2)]])
barr = np.block([[barr],[np.ones((N,1))/4]])
return barr
Gmat = Gmatrix(N)
Kmat = Kmatrix(N,Gmat,Imat)
Barr = Barray(N)
Phi = np.linalg.solve(Kmat,Barr)
def trans(N,Phi):
# 把线性方程组解得的列矢量分割，从左到右排列转化为方矩阵，方便可视化
phi = Phi[0:N]
for i in range(N-1):
j = i+1
phi = np.block([phi,Phi[(j*N):((j+1)*N)]])
return phi
Phi = trans(N,Phi)
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
X,Y = np.meshgrid(x,y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X,Y,Phi,rstride=1,cstride=1,cmap=plt.get_cmap('rainbow'))
ax.contourf(X,Y,Phi,zdir='z',offset=1,cmap='rainbow')
plt.title('Difference methods for Laplace equation',fontsize='12')
plt.show()
作图如下：
通过控制作图参数里的offset，可以调整等高线的位置：
差分法求解的误差主要来自于差分公式对偏微分的近似替代，为了提高求解精度，可以把网格划分得更密，同时计算量也会显著增大。
Reference：
1、马文淦，《计算物理学》
2、李荣华，刘播，《微分方程数值解法》