General Discussion of Schrödinger Equation

single particle equation

 $$i\hbar \frac{\partial}{\partial t }|\psi(\mathbf{x},t)\rangle=H|\psi(\mathbf{x},t)\rangle,H=-\frac{\hbar ^2}{2m}\nabla^2+V(\mathbf{x},t)$$ 




multi-particle

 $$i\hbar \frac{\partial}{\partial t }|\psi(\mathbf{x}_1,\cdots,\mathbf{x}_N,t)\rangle=H|\psi(\mathbf{x}_1,\cdots,\mathbf{x}_N,t)\rangle,H=\sum_{i=1}^N[-\frac{\hbar ^2}{2m}\nabla^2+V(\mathbf{x},t)]+V(\mathbf{x}_1,\cdots,\mathbf{x}_N,t)$$ 




Properties

 $$\frac{\partial}{\partial t}\rho+\nabla\cdot\mathbf{j}=0$$ 


where 

 $$\rho=|\psi(\mathbf{x},t)^2|=\psi(\mathbf{x},t)^*\psi(\mathbf{x},t) \\ \mathbf{j}=\frac{\hbar}{2m\mathrm{i}}[\psi^*(\mathbf{x},t)\nabla\psi(\mathbf{x},t)-\psi(\mathbf{x},t)\nabla\psi^*(\mathbf{x},t)]$$ 



  
proof 

  We assume $V(\mathbf{x},t)$ to be real 

  
 $$i\hbar \frac{\partial}{\partial t }\psi(\mathbf{x},t)=-\frac{\hbar ^2}{2m}\nabla^2\psi(\mathbf{x},t)+V(\mathbf{x},t)\psi(\mathbf{x},t)\\ \Rightarrow -i\hbar \frac{\partial}{\partial t }\psi^*(\mathbf{x},t)=-\frac{\hbar ^2}{2m}\nabla^2\psi^*(\mathbf{x},t)+V(\mathbf{x},t)\psi^*(\mathbf{x},t)\\ \Rightarrow\begin{cases}\psi^*(\mathbf{x},t)i\hbar \dfrac{\partial}{\partial t }\psi(\mathbf{x},t)=\psi^*(\mathbf{x},t)\left(-\dfrac{\hbar ^2}{2m}\nabla^2\psi(\mathbf{x},t)+V(\mathbf{x},t)\psi(\mathbf{x},t)\right)\\  \psi(\mathbf{x},t)\left(-i\hbar \dfrac{\partial}{\partial t }\psi^*(\mathbf{x},t)\right)=\psi(\mathbf{x},t)\left(-\dfrac{\hbar ^2}{2m}\nabla^2\psi^*(\mathbf{x},t)+V(\mathbf{x},t)\psi^*(\mathbf{x},t)\right)\end{cases}$$  

  $(1)-(2)$ 

   $$l.h.s=i\hbar(\psi^*\frac{\partial}{\partial t}\psi+\psi\frac{\partial}{\partial t}\psi^*)=i\hbar \frac{\partial}{\partial t}|\psi(\mathbf{x},t)^2|\\r.h.s=-\dfrac{\hbar ^2}{2m}(\psi^*\nabla^2\psi-\psi\nabla^2\psi^*)\\=-\dfrac{\hbar ^2}{2m}\nabla\cdot(\psi^*\nabla\psi-\psi\nabla\psi^*)\\\Rightarrow i\hbar \frac{\partial}{\partial t}|\psi(\mathbf{x},t)^2|+\frac{\hbar}{2mi}\nabla\cdot(\psi^*\nabla\psi-\psi\nabla\psi^*)\equiv0$$  

  let $\rho=|\psi(\mathbf{x},t)^2|,\mathbf{j}=\dfrac{\hbar}{2m\mathrm{i}}[\psi^*\nabla\psi-\psi\nabla\psi^*]$ 

  thus 

   $$\frac{\partial}{\partial t}\rho+\nabla\cdot\mathbf{j}=0$$ 



Some properties of 1D stationary problem

 $$i\hbar \frac{\partial}{\partial t }|\psi(\mathbf{x},t)\rangle=H|\psi(\mathbf{x},t)\rangle$$  

Assume $|\psi(\mathbf{x},t)\rangle=|\psi(\mathbf{x})\rangle \cdot e^{-iEt/\hbar}$ 

Then 

 $$H\psi=E\psi$$ 



$\left[-\dfrac{\hbar ^2}{2m}\nabla^2+V(\mathbf{x})\right]\psi(\mathbf{x})=E\psi(\mathbf{x})$ 

正问题：$V(\mathbf{x})\rightarrow E,\psi(\mathbf{x})$ 

逆问题：$E \rightarrow V(\mathbf{x})$
both $\psi(\mathbf{x}),\psi^*(\mathbf{x})$ are solutions.
if $\psi(\mathbf{x})$ is not degenerate, it can be chosen to be real.
if $\psi(\mathbf{x})$ is a solution and $V(\mathbf{x})=V(-\mathbf{x})$ then $\psi(-\mathbf{x})$ is another solution.
if $V(\mathbf{x})=V(-\mathbf{x})$, a complete set of solutions.


Infinite Square Potential Well Problem

（上次码好第六次笔记发布之后，不小心被第七次覆盖，简直欲哭无泪。上次只备份到这里，以后有机会再补上吧，sigh……）



Thanks to Prof. Guo Hong

摘要：此文介绍了一种使用MATLAB求解一维定态薛定谔方程的方法。利用充分格式进行离散化，得出相应的矩阵方程，用MATLAB求解本征值和本征函数。此方法简单可靠，可以处理各种时间无关的束缚态问题。所用的程序附于文后。
宏观物体遵循的是牛顿运动方程，给定初始条件以及受力条件，我们就可以求出任意时刻粒子的状态。而在原子尺度上，所有粒子都表现出波的行为，粒子的状态用波函数
来描述，描述微观粒子运动的方程是由薛定谔于1925年提出薛定谔方程。微观粒子的运动与其所处的势场相关，当势场不随时间变化时，称为定态，一维情况下的定态薛定谔方程为
其中， 
 表示粒子所处的势场。由定态波函数可以得出总波函数
1926年， Born提出，波函数模的平方
 代表在位置 
 ，时刻 
 寻找到电子的概率，这就是波函数的条统计解释。可以由遵循微分方程的波函数表示，薛定谔波方程涉及空间坐标和时间。对于一维情况，在时刻 
 时，在 
 和 
 之间的某处找到电子的概率由下式给出
时间无关薛定谔方程(1)是系统的能量本征方程。该特征值方程通常由一组特定的函数 
 和相应的一组常数 
 满足解的条件，被称为是哈密顿算符的本征函数和相应的本征值。当测量处于 
 状态的系统的能量时，结果将始终为 
 。对于束缚态情况，必须有
一维定态薛定谔方程是二阶微分方程，但是，能够解析求解的情况屈指可数，如氢原子，谐振子，无限深势阱等。随着计算机技术的发展，我们可以数值上进行求解。本文利用MATLAB软件，使用矩阵方法求解束缚态的本征值问题。
对于原子系统，以nm和eV的能量测量长度更方便。我们可以使用缩放因子
因此我们可以将等式(1)写成
为了求解上述方程式，首先进行离散化处理。将坐标 
 离散化为 
 个点，用 
 来表示，若 
 ，则有 
 另外，还需要对方程式进行差分格式处理，二阶导数处理如下
因此， 
 的二阶导数矩阵可以写成 
 矩阵
 矩阵大小为 
 而不是 
 ，因为函数的二阶导数无法在终点进行计算，即 
 和 
 。动能矩阵为 
 ，势能矩阵为对角矩阵，即 
 。则我们现在可以将哈密顿矩阵定义为 
 。用于生成哈密顿矩阵的代码是
% 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;
因此，矩阵形式的薛定谔方程是
。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 : num
    if abs(x(cn)) <= x1
        U(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 == 0
    E(n) = e_values(n,n);
    if E(n) > 0
        flag = 1;
    end % if
    n = n + 1;
end  % while
E(n-1) = [];
n = n-2;
% Corresponding Eigenfunctions 1, 2, ... ,n: Normalizing the wavefunction
for cn = 1 : n
    psi(:,cn) = [0; e_funct(:,cn); 0]; 
    area = simpson1d((psi(:,cn) .* psi(:,cn))',xMin,xMax);
    psi(:,cn) = psi(:,cn)/sqrt(area);       % normalize
    prob(:,cn) = psi(:,cn) .* psi(:,cn);
    if psi(5,cn) < 0
        psi(:,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 : n
    fprintf('  %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 : cnmax
  ys(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 < 6
    nMax = n;
else
    nMax = 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:nMax
    subplot(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 on
    plot(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 on
    plot(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) == 1
    sc = 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

 

1.定态薛定谔方程

   薛定谔方程为量子力学中的基本方程，其揭示微观世界物质运动基本规律、是原子物理学中处理一切非相对论问题的有力工具。本文利用基础计算机数值模拟方法对其能级、波函数等性质求解。

   薛定谔方程通常表示：

                                                     

   若求能级与波函数，势函数不含时间，则代入下述方程： 

                                                    

   得到定态不含时的薛定谔方程：

                                                  

   2.一维方势阱中粒子能级、波函数模拟求解

                                                       

      对于如图一维方势阱问题，其势函数为：

                                                        

       带入上述微分方程，其解为：

                                                   

      利用边界条件确定系数A,B,C,D

     利用波函数连续性可得：

 

                                             

      求得：

                                                    

 

     经过推导得到归一化的本征函数：

                                                              

      3.matlab数值计算

     其计算流程图如下：

                                       

N=0;
V0=-20;
W=1.0;%势阱宽度
Emin=V0%势能
Emax=0;
M=51;%离散步长
de=(Emax-Emin)./(M-1);
E=Emin-de;
Ans=[];
for i=1:M-1
    E=E+de;
    r1=sqrt(abs(E-V0));
    r2=sqrt(abs(E));
    phi=(r2./r1)*sin(W*r1)+cos(W*r1);
    phi1=(r2*cos(W*r1)-r1*sin(W*r1));
    C1=0.5*(exp(-W*r2))*(phi+(phi1./r2));
    D1=0.5*(exp(W*r2))*(phi-(phi1./r2));
    N=floor((W*r1)./3.1415926);
    if (N-floor(0.5*N)*2~=0&&phi>0)
        N=N+1;
    end
    if (N-floor(0.5*N)*2==0&&phi<0)
        N=N+1;
    end
    if C1*D1<0
        if ((log(-D1./C1))./(2*r2))>W
            N=N+1;
            Ans(i,1)=E;
            Ans(i,2)=N;
        else
            Ans(i,1)=E;
            Ans(i,2)=N;
        end
    else
        Ans(i,1)=E;
        Ans(i,2)=N;
    end
end
count=2;
n=[0];%N值量子数
E=[Emin];
for i=1:size(Ans(:,2),1)
    if Ans(i,2)~=n(count-1)
        n(count)=Ans(i,2);
        E(count)=Ans(i,1);
        count=count+1;
    end
end
plot(n,E,'ro');
%ylim([-3,3]);
title('E与n对应值');
xlabel('n')
ylabel('E')

计算结果如下：

                           

                                                                            表 变量E与N对应值

1


			
E
			
			

			
-20
			
			

			
…
			
			

			
-15.6
			
			

			
-15.2
			
			

			
-14.8
			
			

			
…
			
			

			
-4.0
			
			

			
-3.6
			
			

			
…
			
			

			
-0.4
			
		

			
N
			
			

			
0
			
			

			
…
			
			

			
0
			
			

			
1
			
			

			
1
			
			

			
…
			
			

			
1
			
			

			
2
			
			

			
…
			
			

			
2
			
		
 

波函数绘图

 

 

L1=5;
x=0:0.01:L1;%势阱宽度
z1=(sqrt(2./L1))*sin((1*pi*x)./L1);
z2=(sqrt(2./L1))*sin((2*pi*x)./L1);
z3=(sqrt(2./L1))*sin((3*pi*x)./L1);
subplot(3,2,1);
plot(x,z1.^2);
%hold onlegend('n=1','n=2','n=3')
legend('n=1')
title('波函数');
xlabel('势阱宽度')
ylabel('概率密度')

%hold onlegend('n=1','n=2','n=3')
%%%%%%%%%%%%
subplot(3,2,2);
plot(x,z1);

legend('n=1')
title('波函数');
xlabel('势阱宽度')
ylabel('|波函数|')

%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,2,3);
plot(x,z2.^2);
legend('n=2')
title('波函数');
xlabel('势阱宽度')
ylabel('概率密度')
subplot(3,2,4);
plot(x,z2.^2);
legend('n=2')
title('波函数');
xlabel('势阱宽度')
ylabel('概率密度')
%%%%%%%%%%%%%
%hold on
subplot(3,2,5);
plot(x,z3.^2);
legend('n=3')
title('波函数');
xlabel('势阱宽度')
ylabel('概率密度')
subplot(3,2,6);
plot(x,z3.^2);
legend('n=3')
title('波函数');
xlabel('势阱宽度')
ylabel('|波函数')
%set(gca,'ytick',[],'ycolor','w')


 

 

			


伯努利方程 

Bernoulli equation
 
流体宏观运动机械能守恒原理的数学表达式。1738年瑞士数学家D.伯努利在《水动力学──关于流体中力和运动的说明》中提出了这一方程。它可由理想流体运动方程（即欧拉方程）在定态流动条件下沿流线积分得出；也可由热力学第一定律导出。它是一维流动问题中的一个主要关系式，在分析不可压缩流体的定态流动时十分重要，常用于确定流动过程中速度和压力之间的相互关系。

 
方程的形式  
对于不可压缩的理想流体，密度不随压力而变化，可得： 
 


式中Z为距离基准面的高度;p为静压力;u为流体速度;ρ为流体密度;g为重力加速度。方程中的每一项均为单位质量流体所具有的机械能，其单位为N·m/kg，式中左侧三项，依次称为位能项、静压能项和动能项。方程表明三种能量可以相互转换，但总和不变。当流体在水平管道中流动时Z不变，上式可简化为：


 
此式表述了流速与压力之间的关系:流速大处压力小,流速小处压力大。 
 
对于单位重量流体,取管道的1、2两截面为基准,则方程的形式成为：


 
式中每一项均为单位重量流体的能量，具有长度的因次，三项依次称为位头、静压头和动压头（速度头）。 
 
对于可压缩理想流体，密度随压力而变化。若这一变化是可逆等温过程，则方程可写成下式： 
 


若为可逆绝热过程，方程可写为： 
 


式中γ为定压比热容cp和定容比热容cV之比，即比热容比，也称为绝热指数。 
 
对于粘性流体，流动截面上存在着速度分布，如用平均流速ū表达动能项，应对其乘以动能校正系数

α。此外,还需考虑因粘性引起的流动阻力,即造成单位质量流体的机械能损失hf，若在流体流动过程中，单位质量流体又接受了流体输送机械所做的功W，在这些条件下,若取处于均匀流段的两截面1和2为基准，则方程可扩充为： 
 


α

值可由速度分布计算而得, 
流体在圆管内作层流流动时 
α=2；作湍流流动时，

α≈1.06。 
 
方程的应用  
伯努利方程阐明的位能、动能、静压能相互转换的原理，可用来分析计算一些实际问题，例如： 
 
①计算流体从小孔流出的流速
设在容器中盛有液体，液面维持不变,距液面下h处的容器壁面上开有一小孔，液体在重力作用下自小孔流出。据伯努利方程可以计算出液体由小孔流出时的平均流速为： 
 


式中Cd为孔流系数，其值由实验确定，约为0.61～0.62；g为重力加速度。由上述速度及已知的小孔面积,可算出通过小孔的流量；或由这一关系，计算确定达到一定流量所必须维持的液面高度。若气体在一定压力差作用下由容器壁上的小孔流出，当速度不过大时，可视为不可压缩流体，其流量也可以利用伯努利方程来估计。 
 


②毕托管 设均匀气流以等速uo绕过某物体流动，气流受阻后在物体前缘（A处）停滞，形成驻点（图1），该点处的压力称为驻点压力pA。若未受扰动的某点O压力为po，由伯努利方程可得 
 


测出pA与po的差值,即可算出流速uo。据此原理计设的测速装置，称测速器,又称毕托管。毕托管（图2）由一个圆头的双层套管组成，在圆头中心处开有与内套管相连的小孔,内套管与测压计的一头联接,以测定驻点压力pA；在外套管侧表面一定距离处，沿周向均匀地开一排与管壁垂直的静压孔，外套管与测压计的另一头相联，以测定压力po。根据测得的压力差h,可计算测点处的流速。 
 


③文丘里管 又称文氏管(图3)，是一种先收缩而后逐渐扩大的管道。由于截面积有变化，流速改变，根据伯努利方程，压力也随之改变。量出管前与喉管处的压力差，即可推算流量。用于测量流量的文丘里管，称文丘里流量计。又由于文丘里管喉部形成高速气流，会产生负压而抽吸液体，使气液密切接触，用于完成气体的洗涤、冷却、吸收和反应等操作。用于这类操作的文丘里管称为文丘里洗涤器。


伯努利方程 
Bernoulli
equation
 
流体在忽略粘性损失的流动中，流线上任意两点的压力势能、动能与位势能之和保持不变。这个理论是由瑞士数学家丹尼尔第一·伯努利在1738年提出的，当时被称为伯努利原理。后人又将重力场中欧拉方程在定常流动时沿流线的积分称为伯努利积分，将重力场中无粘性流体定常绝热流动的能量方程称为伯努利定理。这些统称为伯努利方程，是流体动力学基本方程之一。

 
重力场中无粘性不可压缩流体定常流动的伯努利方程有以下3种形式

 (1)


  (2)


  (3)


式中下角标"1"、"2"分别表示流线上任意两点的位置；式(1)中、、gz和E分别视为单位质量流体的压力势能、动能、位势能和总能量；式(2)常用在水力工程方面。式中、、z 
和H 
分别称为压力水头、速度水头、位置水头和总水头;式(3)经常用于流体测量方面，它忽略了重力项的影响。式中p、和p0分别称为静压、动压和总压。以上
3式分别表示流线上任意两点的总能量、总水头和总压相同。

 
图中的水箱小孔出流，可作为伯努利方程的应用实例。p1=p2=pα，pα为大气压力，因为流速，由式(1)可得小孔出流速度。其中h为小孔到液面的高度。





对于重力场中无粘性气体的定常绝热流动，伯努利方程有以下形式 
  
式中h1、h2分别为流线上任意1、2两点气体的比焓（见焓）；h 0为气体的总焓。上式表示流线上任意两点见煮沸法,CaCO3,MgCO3,药剂软化法,离子交换法,电渗析和超滤技术焓相等。对于重力场中粘性气体定常绝热管内流动，上式也适用，这时下角标1、2表示任意两横截面的位置，v1、v2、h1、h2、h0均取对应管道横截面上的平均值。

 
							
		
转载于:https://www.cnblogs.com/crid/archive/2012/09/21/5172525.html