PHYS 410辅导 、讲解 MATLAB程序设计
PHYS 410: Computational Physics Fall 2024
Project 2—The Time Dependent Schro¨dinger Equation
Due: Monday, December 2, 11:59 PM
1. Problem 1
1.1 Introduction
In this part of the project you will solve the one-dimensional time-dependent Schr¨odinger equation using the method
discussed in the lecture, which will be summarized here.
The continuum equation is (after non-dimensionalization)
iψ(x,t)t = −ψxx + V (x,t)ψ , (1)
where the wavefunction, ψ(x,t), is complex. The equation is to be solved on the domain
0 ≤ x ≤ 1 , 0 ≤ t ≤ tmax,
subject to initial and boundary conditions
ψ(x, 0) = ψ0(x), (2)
ψ(0,t) = ψ(1,t) = 0 . (3)
The initial data function ψ0(x) can be complex in general, but in the cases we will consider will be real. We will also
restrict attention to time-independent potentials, but carry the explicit dependence in the following to emphasize
that, in principle, it is no more difficult to deal with a time-dependent potential than a time-independent one.
We note that ψ can be expressed in terms of its real and imaginary parts, ψRe and ψIm, respectively as
ψ = ψRe + iψIm
A useful diagnostic quantity is the “running integral”, P(x,t), of the probability density, ρ = |ψ|
2 = ψψ
∗
:
P(x,t) =
Z x
0
ψ(x˜,t)ψ
⋆
(x˜,t)dx˜ (4)
If the wavefunction is properly normalized, then we will have
P(1,t) = 1
but even if it is not so normalized (and in our applications there will be no need to ensure normalization), we should
have
P(1,t) = conserved to level of solution error
The quantity
√ρ = |ψ| is also a useful diagnostic. However, as you develop your code to solve the Schr¨odinger
equation you should be prepared to plot the real and imaginary parts of ψ as well.
A family of exact solutions of (1) is
ψ(x,t) = e
−im2π
2
t
sin(mπx) (5)
where m is a positive integer.
1.2 Discretization
We discretize the continuum domain by introducing the discretization level, l, and the ratio of temporal to spatial
mesh spacings
λ =
∆t
∆x
1Then
nx = 2
l + 1
∆x = 2
−l
∆t = λ∆x
nt = round (tmax/∆t) + 1
We apply the Crank-Nicolson discretization approach to (1)–(3) yielding
i
ψ
n+1
j
− ψ
n
j
∆t
= −
1
2
ψ
n+1
j+1
− 2ψ
n+1
j
+ ψ
n+1
j−1
∆x
2 +
ψ
n
j+1
− 2ψ
n
j
+ ψ
n
j−1
∆x
2
+
1
2
V
n+ 1
2
j
ψ
n+1
j
+ ψ
n
j
j = 2, 3, . . . , nx − 1 , n = 1, 2, . . . nt − 1 (6)
ψ
n+1
1
= ψ
n+1
nx
= 0 , n = 1, 2, . . . nt − 1 (7)
ψ
1
j
= ψ0(xj ), j = 1, 2, . . . nx (8)
1.3 Implementation
Implement your solution of (6)–(8) as a MATLAB function with the following header and arguments:
function [x t psi psire psiim psimod prob v] = ...
sch_1d_cn(tmax, level, lambda, idtype, idpar, vtype, vpar)
% Inputs
%
% tmax: Maximum integration time
% level: Discretization level
% lambda: dt/dx
% idtype: Selects initial data type
% idpar: Vector of initial data parameters
% vtype: Selects potential type
% vpar: Vector of potential parameters
%
% Outputs
%
% x: Vector of x coordinates [nx]
% t: Vector of t coordinates [nt]
% psi: Array of computed psi values [nt x nx]
% psire Array of computed psi_re values [nt x nx]
% psiim Array of computed psi_im values [nt x nx]
% psimod Array of computed sqrt(psi psi*) values [nt x nx]
% prob Array of computed running integral values [nt x nx]
% v Array of potential values [nx]
The input parameters idtype and vtype are integers that select which initial data type and potential type, respectively,
are to be used. Dependent on the type, elements of the associated parameter vector will be used to define the
initial data or potential. Specifically, you are to implement options as follows:
Initial data types
• idtype == 0: Exact family (5)
ψ(x, 0) = sin(mπx)
– idpar(1): m
2• idtype == 1: Boosted Gaussian:
ψ(x, 0) = e
ipx
e
−((x−x0)/δ)
2
– idpar(1): x0
– idpar(2): δ
– idpar(3): p
Potential types
• vtype == 0: No potential.
V (x) = 0
• vtype == 1: rectangular barrier or well.
V (x) =
0 for x < xmin
Vc for xmin ≤ x ≤ xmax
0 for x > xmax
– vpar(1): xmin
– vpar(2): xmax
– vpar(3): Vc
Notes
• The terminology “boosted Gaussian” comes from the fact that by pre-multiplying the (real-valued) Gaussian
profile by e
ipx
, we give the wave packet some momentum in the direction of p (which can have either sign).
• The constant Vc is positive for a barrier, negative for a well.
• Don’t worry about the fact that, strictly speaking, the boosted Gaussian profile is incompatible with the
boundary conditions. In practice we will try to center the Gaussian sufficiently far from the boundaries that
the incompatibility is lost in the truncation error. You should always set the wave function to 0 at x = 0 and
x = 1, including at the initial time.
Coding
• Write equations (6)–(7) as a complex (as in complex number) tridiagonal system for the advanced unknowns
ψ
n+1
j
. In your code, set up the tridiagonal system using spdiags as discussed in class and as illustrated in the
code diff 1d imp.m from Tutorial 6. and solve it using left division. Note that MATLAB has native support for
complex numbers (the variables i and j are both initialized to the unit pure imaginary number, i) and that
you should implement your solution directly using complex arithmetic. As observed above, there is no need to
worry about normalization of the wave function.
• Once you’ve set up the tridiagonal matrix, say A, using spdiags, you can view the corresponding full matrix
using the MATLAB command full(A). This can be useful while debugging.
• Note that the MATLAB abs(z) command when applied to a complex number, z, returns its modulus, |z|.
• The integral in (4) can be computed to O(h
2
) using the trapezoidal formula. Recall that if we are given n
approximate values fi at values of x, xi
, then the trapezoidal approximation is given by
Z xn
x1
f(x)dx ≈
1
2
nX−1
i=1
(fi + fi+1)(xi+1 − xi)
1.4 Convergence testing
Define a level l solution computed using sch 1d cn by ψ
l
. Note that ψ
l
is a function of both the discrete time and
space coordinates. Denote by dψ
l
the quantity defined by
dψ
l = ψ
l+1 − ψ
l
3where it is to be understood that the data defined by the grid function (array) ψ
l+1
is 2:1 coarsened in both the
time and space dimensions so that it has the same size as ψ
l
. Then one way we can convergence test sch 1d cn is to
compute
kdψ
l
k2(t
n
) (9)
where
k · k2
denotes the l-2 norm (RMS value) that has been defined before; namely, for any length-m vector v
kvk2 =
sPm
j=1
|v
j
|
2
m
Observe that for complex numbers, |v
j
| is the modulus of the number. Note that (9) involves taking spatial norms
of the pairwise subtraction of grid functions at two different levels. This results in a function of the discrete time,
t
n, on the level-l grid.
Following the development we have seen for solutions of other finite difference equations, we note that since our FDA
is O(h
2
) (where ∆x = h and ∆t = λh), we expect the solution to be O(h
2
) accurate, with ψ
l admitting an expansion
of the form
ψ
l
(x,t) = ψ(x,t) + h
2
l
e2(x,t) + O(h
4
l
) (10)
where e2(x,t) is some error function. From (10) we can deduce that if we graph rescaled values of kdψ
lk2 on a single
plot, then convergence is signalled by near-coincidence of the curves, with better agreement as we go to higher values
of l. In particular, for a test with levels
l = lmin, lmin + 1, . . . , lmax
we should plot
kdψ
lmin k2 , 4kdψ
lmin+1
k2 , 4
2
kdψ
lmin+2
k2 , . . . , 4
lmax−lmin−1
kdψ
lmax−1
k2
Additionally, for the case when idtype = 0, so that we know the exact solution, we can compute the actual solution
errors. Specifically, we can perform precisely the same type of convergence test just described, but where ψ
l+1
is
replaced with ψexact. Thus we define
kE(ψ
l
)k2(t
n
) = kψexact − ψ
l
k2(t
n
)
and use exactly the same plotting strategy for kE(ψ
l
)k2 as we do for kdψ
lk2.
Convergence tests to perform
Ensure that at least one of the following convergence tests, including the plotting, can be performed by executing a
script ctest 1d.
1. idtype = 0, vtype = 0
Other parameters
• idpar = [3]
• tmax = 0.25
• lambda = 0.1
• lmin = 6
• lmax= 9
Perform the 4-level test and make convergence plots as described above for both kdψ
lk2 and kE(ψ
l
)k2. Include
the plots in your report.
2. idtype = 1, vtype = 0
Other parameters
• idpar = [0.50 0.075 0.0]
• tmax = 0.01
4• lambda = 0.01
• lmin= 6
• lmax= 9
Perform the 4-level test and make a convergence plot as described above for kdψ
lk2. Include the plot in your
report.
1.5 Numerical Experiments
Consider the discrete running integral of the probability density:
P
n
j
= P(xj ,t
n
), j = 1, 2, . . . nx , n = 1, 2, . . . nt
Define the temporal average, P¯
j
, of the above quantity:
P¯
j
=
Pnt
n=1 P
n
j
nt
and note that the MATLAB command mean can be used to compute this. We will also want to ensure that P¯
j
is
properly normalized so that P¯
nx = 1. We can do this as follows:
P¯
j
:= P¯
j
/P¯
nx
, j = 1, 2, . . . nx
In the following we will assume that P¯
j
has been properly normalized. Given two values of x, x1 and x2, satisfying
x2 > x1, we can interpret the quantity
P¯(x2) − P¯(x1)
as the fraction of time our quantum particle spends in the interval x1 ≤ x ≤ x2. Here the notation P¯(x) is to
be interpreted as P¯(x) = P¯(xj ) = P¯
j
where xj is the nearest grid point to x. Now, for a free particle—i.e. for
V = 0—and for sufficiently long times, we expect that
P¯(x2) − P¯(x1) → x2 − x1
That is, the fraction of time the particle spends in the interval is given simply by the width of the interval (this
direct equality is due to the fact that we are solving the Schr¨odinger equation on the unit interval, 0 ≤ x ≤ 1).
For the general case of a non-zero potential, we can then define the excess fractional probability that the particle
spends in a given spatial interval as
F¯
e(x1, x2) =
P¯(x2) − P¯(x1)
x2 − x1
For the experiments described below, this quantity will span orders of magnitude so, particularly for the purposes of
plotting, it will be convenient to compute its (natural) logarithm (recall that MATLAB uses log for the natural log).
ln F¯
e(x1, x2) = ln
P¯(x2) − P¯(x1)
x2 − x1
Also observe that although we have called F¯
e(x1, x2) the excess fractional probability, it can in fact satisfy F¯
e(x1, x2) <
1, indicating that the particle is spending less time in the specified interval than a free particle would. Indeed, in the
experiments that follow, you should find that F¯
e(x1, x2) < 1 is the rule rather than the exception.
51.5.1 Experiment 1: Barrier Survey
In this investigation the particle will start to the left of a barrier whose height, V0, is the control parameter for the
experiment. You will then determine the dependence of ln(F¯
e(x1, x2)) on ln(V0) where x1 and x2 span the region to
the right of the barrier.
Parameters
• tmax = 0.10
• level = 9
• lambda = 0.01
• idtype = 1 (boosted Gaussian)
• idpar = [0.40, 0.075, 20.0]
• vtype = 1 (rectangular barrier)
• vpar = [0.6, 0.8, VARIABLE > 0]
• x1 = 0.8
• x2 = 1.0
Write a script called barrier survey that uses sch 1d cn to compute F¯
e(0.8, 1.0) for 251 uniformly spaced values
of ln(V0) ranging from -2 to 5. The script is to make a plot of ln(F¯
e(0.8, 1.0)) versus ln(V0). Include the plot in your
writeup and comment on what you can deduce from it. Since it will take some time for the 251 runs to complete, it
only makes sense to debug your script using fewer values of ln(V0).
1.5.2 Experiment 2: Well Survey
The second experiment is very much like the first, except that we now consider scattering of a particle off a potential
well. Once again you will perform a survey to investigate the dependence of ln(F¯
e(x1, x2)) on ln(V0) where x1 and
x2 span the location of the well.
Parameters
• tmax = 0.10
• level = 9
• lambda = 0.01
• idtype = 1 (boosted Gaussian, but note that p = 0)
• idpar = [0.40, 0.075, 0.0]
• vtype = 1 (rectangular well)
• vpar = [0.6, 0.8, VARIABLE < 0]
• x1 = 0.6
• x2 = 0.8
Write a script called well survey that uses sch 1d cn to compute F¯
e(0.6, 0.8) for 251 uniformly spaced values of
ln(|V0|) ranging from 2 to 10. Note that V0 is strictly less than 0. The script is to make a plot of ln(F¯
e(0.6, 0.8))
versus ln(|V0|). Include the plot in your writeup and discuss what you can conclude from it.
62. Problem 2
2.1 Introduction
In this problem you will solve the two-dimensional Schr¨odinger equation using the ADI technique discussed in class
for the case of the diffusion equation. Indeed, as mentioned in lecture, and should you have time, I suggest that you
consider first solving the 2d diffusion equation using ADI which will give you a good base from which to tackle the
Schr¨odinger equation.
The non-dimensionalized continuum equation is
iψ(x, y, t)t = −(ψxx + ψyy) + V (x, y)ψ
or, multiplying through by −i
ψt = i(ψxx + ψyy) − iV (x, y)ψ (11)
In this last form, the similarity to a diffusion equation with an imaginary diffusion constant (and a source term) is
apparent. Equation (11) is to be solved on the domain
0 ≤ x ≤ 1 , 0 ≤ y ≤ 1 , 0 ≤ t ≤ tmax
subject to initial and boundary conditions
ψ(x, y, 0) = ψ0(x, y) (12)
ψ(0, y, t) = ψ(1, y, t) = ψ(x, 0, t) = ψ(x, 1, t) = 0 (13)
A family of exact solutions is given by
ψ(x, y, t) = e
−i(m2
x+m2
y
)π
2
t
sin(mxπx) sin(myπy) (14)
where mx and my are positive integers.
2.2 Discretization and the ADI scheme
As for the 1d case, the continuum domain is discretized by introducing the discretization level, l, and the ratio of
temporal to spatial mesh spacings
λ =
∆t
∆x
=
∆t
∆y
Then
nx = ny = 2l + 1
∆x = ∆y = 2−l
∆t = λ∆x
nt = round (tmax/∆t) + 1
Defining the difference operators ∂
h
xx and ∂
h
yy by
∂
h
xxu
n
i,j ≡
u
n
i+1,j − 2u
n
i,j + u
n
i−1,j
∆x
2
∂
h
yyu
n
i,j ≡
u
n
i,j+1 − 2u
n
i,j + u
n
i,j−1
∆y
2
the following is an ADI discretization of (11):
1 − i
∆t
2
∂
h
xx
ψ
n+ 1
2
i,j =
1 + i
∆t
2
∂
h
xx 1 + i
∆t
2
∂
h
yy − i
∆t
2
V
i,j
ψ
n
i,j ,
i = 2, 3, . . . , nx − 1 , j = 2, 3, . . . , ny − 1 , n = 1, 2, . . . , nt − 1 (15)
7
1 − i
∆t
2
∂
h
yy + i
∆t
2
V
i,j
ψ
n+1
i,j = ψ
n+ 1
2
i,j ,
i = 2, 3, . . . , nx − 1 , j = 2, 3, . . . , ny − 1 , n = 1, 2, . . . , nt − 1 (16)
Equations (15) and (16) are to be supplemented with the initial conditions
ψ
1
i,j = ψ0(xi
, yj ), (17)
and the boundary conditions
ψ
n
1,j = ψ
n
nx,j = ψ
n
i,1 = ψ
n
i,ny = 0 . (18)
2.3 Implementation
Implement your solution of (15)–(18) as a MATLAB function, sch 2d adi with the following header and arguments:
function [x y t psi psire psiim psimod v] = ...
sch_2d_adi(tmax, level, lambda, idtype, idpar, vtype, vpar)
% Inputs
%
% tmax: Maximum integration time
% level: Discretization level
% lambda: dt/dx
% idtype: Selects initial data type
% idpar: Vector of initial data parameters
% vtype: Selects potential type
% vpar: Vector of potential parameters
%
% Outputs
%
% x: Vector of x coordinates [nx]
% y: Vector of y coordinates [ny]
% t: Vector of t coordinates [nt]
% psi: Array of computed psi values [nt x nx x ny]
% psire Array of computed psi_re values [nt x nx x ny]
% psiim Array of computed psi_im values [nt x nx x ny]
% psimod Array of computed sqrt(psi psi*) values [nt x nx x ny]
% v Array of potential values [nx x ny]
As before, the integer arguments idtype and vtype encode the initial data and potential types, respectively, with
idpar and vpar defining the parameters for the various options. In this case you are to implement options as follows:
Initial data types
• idtype == 0: Exact family (14)
ψ(x, y, 0) = sin(mxπx) sin(myπy)
– idpar(1): mx
– idpar(2): my
• idtype == 1: Boosted Gaussian
ψ(x, y, 0) = e
ipxx
e
ipyy
e
−((x−x0)
2/δ2
x+(y−y0)
2/δ2
y
)
– idpar(1): x0
– idpar(2): y0
– idpar(3): δx
– idpar(4): δy
8– idpar(5): px
– idpar(6): py
As previously, don’t worry about the fact that the Gaussian data is strictly speaking incompatible with the
boundary conditions; just be sure to always impose the correct boundary conditions.
Potential types
• vtype == 0: No potential.
V (x, y) = 0
• vtype == 1: Rectangular barrier or well.
V (x, y) =
Vc for (xmin ≤ x ≤ xmax) and (ymin ≤ y ≤ ymax)
0 otherwise
– vpar(1): xmin
– vpar(2): xmax
– vpar(3): ymin
– vpar(4): ymax
– vpar(5): Vc
• vtype == 2: Double slit. Let j
′ = (ny − 1)/4 + 1. Then
Vi,j′ = Vi,j′+1 = 0 for [(x1 ≤ xi) and (xi ≤ x2)] or [(x3 ≤ xi) and (xi ≤ x4)]
Vi,j′ = Vi,j′+1 = Vc otherwise
Vi,j = 0 for j 6= (j
′
or j
′ + 1)
That is, V is only non-zero for y-locations given by j = j
′ or j = j
′ + 1 and for x-positions not coincident with
one of the slits. This thus simulates a thin (two mesh points wide) plate at a fixed y-position, with adjustable
slit openings, which span (x1, x2) and (x3, x4).
– vpar(1): x1
– vpar(2): x2
– vpar(3): x3
– vpar(4): x4
– vpar(5): Vc
Coding
• IMPORTANT!! Beware that the transpose operator ’ computes the conjugate transpose (i.e. it will complexconjugate
any complex numbers it encounters), so probably will not be what you want. Use the non-conjugating
operator .’ instead.
2.4 Convergence testing
Convergence test your code in the same manner that you did for the 1d case, taking into account the extra dimension.
In particular, when performing a convergence test, compute the two-level deviation norms
kdψl
k2(t
n
) (19)
as well as the deviations from the exact solution (when using initial data corresponding to an exact solution)
kE(ψ
l
)k2(t
n
) = kψexact − ψ
l
k2(t
n
)
Note that for the 2d case, the spatial norm k · k2 involves a sum over both spatial dimensions, i.e. treat any 2d grid
function as a vector of length nx × ny.
9Convergence test to perform
• idtype = 0, vtype = 0
Other parameters
– idpar = [2, 3]
– tmax = 0.05
– lambda = 0.05
– lmin = 6
– lmax = 9
Write a script called ctest 2d that performs the 4-level test and makes convergence plots as described above for
both kdψlk2 and kE(ψ
l
)k2. Include the plots in your report.
2.5 Numerical experiments
In the true spirit of projects, you’ll be mostly on your own here. At a minimum, make AVI movies of the following
scenarios and include them in your submission. In all cases you should use Gaussian initial data, with or without a
boost (in either and/or both directions) as you find convenient.
1. Scattering off a rectangular barrier (vtype = 1).
2. Scattering off a rectangular well (vtype = 1).
3. Scattering through a double slit (vtype = 2). Try to simulate discernible particle self-interference on the far
side of the slit.
I recommend that you consider using filled contour plots generated with the MATLAB contourf command as the
basis for making your movies. A useful attribute for contourf is the colormap which can be manipulated using the
colormap command.
Final caution
Note that because your routine will be storing the entire time evolution of the solution (as well as several additional
functions), memory requirements can get extreme for long run times. MATLAB will complain if you try to allocate an
array that is “too large” but you may find your machine gets very sluggish before that point, so beware. At level 9
and below you shouldn’t have too many difficulties, again provided that you don’t try integrating to very long times.
Report any undue difficulties to me as usual.
10