|
Example - Schrödinger Equation
* The "\" symbol has been used at the end of lines to indicate that the line was too long for the display and has been
wrapped to the line below in order to fit within the width of the web page. These lines would need to be joined in the normal
program editor.
I. Introduction
This program demonstrates the use of matrix methods to obtain bound-state solutions to the famous Schrödinger wave equation in quantum mechanics.
II. Mathematics Involved
A 1-dimensional bound state is a time-independent solution ψi to the one-dimensional Schrödinger wave equation  (for suitable units), where V is the potential, Ei are the energy values, which are bounded, vanishes at ± Error! Objects cannot be created from editing field codes. and provides a discrete spectrum of real energy values.
The technique used in the program to solve this is the replacement of this (continuous) differential equation by a (discrete) matrix equation. x was replaced with a set of 2N points, separated by a distance ε ( Nε=1 for the range was used), and ψ(x) was replaced by a vector en. This results in a discretised Schrödinger equation  . M is a 2N by 2N tri-diagonal matrix with leading diagonal entries  , and –0.5 in the diagonals either side of the leading diagonal. The components in the eigenvector ei corresponding to each energy value can be plotted against the components in xi to give a graph of the eigenfunction.
III. How the Program Works
N =15 was used since for this degree of accuracy the eigenvector calculations are still reasonably fast. Increasing N could obtain more precise results, and this could be easily done.
1) The MatSetUp() Procedure
In the procedure MatSetUp() , the elements of the matrix M are assigned their values:
SUB MatSetUp()
FOR i%% = 0 TO 2 * N%% - 1
M.Elements(i%% + 1,i%% + 1) = 1 + epsilon% ^ 2 * V%(x.Elements\
(i%% + 1))
IF i%% > 0 THEN M.Elements(i%%,i%% + 1) = - 0.5
IF i%% < 2 * N%% - 1 THEN M.Elements(i%% + 2,i%% + 1) = - 0.5
NEXT i%%
END SUB
2) The DispEig() Procedure
When the Show button is clicked to display the jth eigensolution, the procedure DispEig() is called which itself calls Plot() :
SUB DispEig(j%%)
'display the jth eigensolution
MOUSE OFF
Superbase.StatusText = "Calculating..."
DIM yscale%
e = M.Eigenvectors(j%%,1)
Superbase.StatusText = "Plotting graph..."
CALL GraphClear()
CALL MaxMin(yscale%)
FOR i%% = 2 TO 2 * N%%
CALL Plot(x.Elements(i%% - 1),e.Elements(i%% - 1),x.Elements(i%%)\
,e.Elements(i%%),yscale%)
NEXT i%%
f.Refresh()
Superbase.StatusText = "Click on an option"
MOUSE ON
END SUB
SUB Plot(xold%,eold%,x%,e%,yscale%)
DIM x1%%,y1%%,x2%%,y2%%
SET fc = f.Add("line" + LTRIM$ ( STR$ (xold%,"99999.00")),"LINE")
fc.ForeColor = 2
x1%% = xs%% + 245 + 225 * xold%:y1%% = ys%% + 195 - yscale% * eold%
x2%% = xs%% + 245 + 225 * x%:y2%% = ys%% + 195 - yscale% * e%
fc.Move(x1%%,y1%%,x2%%,y2%%)
fc.Visible = - 1
END SUB
3) The MaxMin() Procedure
The jth eigenvector is stored in e, then the scale of the plot is calculated in the MaxMin() procedure:
SUB MaxMin(yscale%)
'finds the maximum and minimum values of the eigenfunctions
'to ensure that the function stays on the graph!
DIM i%%
maxval% = Superbase.Maximum("e.Elements(i%%)","i%%",1,2 * N%%)
minval% = Superbase.Minimum("e.Elements(i%%)","i%%",1,2 * N%%)
IF ABS (maxval%) > ABS (minval%) THEN
yscale% = 175 / ABS (maxval%)
ELSE
yscale% = 175 / ABS (minval%)
END IF
END SUB
SUB GraphClear()
FOR xold% = - N%% * epsilon% TO N%% * epsilon% STEP epsilon%
IF f.Exists("line" + LTRIM$ ( STR$ (xold%,"99999.00"))) = - 1 THEN
f("line" + LTRIM$ ( STR$ (xold%,"99999.00"))).Delete()
END IF
NEXT xold%
IF f.Exists("line0.00") = - 1 THEN f("line0.00").Delete()
f.Refresh()
END SUB
The plot is then displayed on screen:
The potential (which is set at first to  , the potential for a harmonic oscillator) can be altered in a dialog box.
|
