Home
|
Course Outline
|
Lectures
|
Homework
|
Files
The program works as follows:
For more information on solving two-point boundary value problems like
this one, see
Chapter 17 of the "Numerical Recipes" textbook.
Here is the code which generates this applet:
// Eigen.java
import comphys.graphics.*;
import comphys.Easy;
public class Eigen extends Animation {
// define potential function
double b = 0.1; // anharmonic coefficient
double V (double x) {
return 0.5 * x * x + b * x * x * x;
}
double E0 = 0.1; // energy guess
double E = E0; // current energy
double dE0 = 0.1; // step for eigenvalue search
double dE = dE0; // current step size
double accuracy = 0.01; // accuracy in energy search
double xLeft = -2; // left integration end point
double xRight = 2; // right integration end point
double xMatch = 0.2; // match point
double dx = 0.02; // integration step size
double[] phiLeft; // wave function in left region
double[] phiRight; // wave function in right region
double[] dPhiLeft; // d phi / dx in left region
double[] dPhiRight; // d phi / dx in right region
int nLeft; // number of points in left region
int nRight; // number of points in right region
int iLeft, iRight; // current indices in integration
double maxPhi; // maximum value of phi
double minPhi; // minimum value of phi
boolean oddNodes; // phi has odd number of nodes
boolean phiMatched; // phiLeft = phiRight at match point
double misMatch; // discrepancy in derivatives at match
double previousMisMatch; // discrepancy at previous energy
boolean eigenFound; // we are done
void initializePhi () {
nLeft = (int) Math.round((xMatch - xLeft) / dx);
if (phiLeft == null || phiLeft.length < nLeft + 1) {
phiLeft = new double[nLeft + 1];
dPhiLeft = new double[nLeft + 1];
}
phiLeft[0] = 0;
dPhiLeft[0] = 1e-10;
for (int i = 1; i <= nLeft; i++)
phiLeft[i] = dPhiLeft[i] = 0;
iLeft = 0;
nRight = (int) Math.round((xRight - xMatch) / dx);
if (phiRight == null || phiRight.length < nRight + 1) {
phiRight = new double[nRight + 1];
dPhiRight = new double[nRight + 1];
}
phiRight[0] = 0;
dPhiRight[0] = -1e-10;
if (oddNodes)
dPhiRight[0] = - dPhiRight[0];
for (int i = 1; i <= nLeft; i++)
phiLeft[i] = dPhiLeft[i] = 0;
iRight = 0;
maxPhi = 0;
minPhi = 0;
oddNodes = phiMatched = false;
}
void integrationStep () {
double phi = 0;
if (iLeft < nLeft) { // left integration not done
double dx = (xMatch - xLeft) / nLeft;
double x = xLeft + iLeft * dx;
phi = phiLeft[iLeft];
double dPhi = dPhiLeft[iLeft];
double d2Phi = 2 * (V(x) - E) * phi;
// use Euler-Cromer algorithm
dPhi += d2Phi * dx;
phi += dPhi * dx;
++iLeft;
phiLeft[iLeft] = phi;
dPhiLeft[iLeft] = dPhi;
if (phiLeft[iLeft] * phiLeft[iLeft - 1] < 0)
oddNodes = !oddNodes;
} else if (iRight < nRight) {
double dx = (xMatch - xRight) / nRight;
double x = xRight + iRight * dx;
phi = phiRight[iRight];
double dPhi = dPhiRight[iRight];
double d2Phi = 2 * (V(x) - E) * phi;
dPhi += d2Phi * dx;
phi += dPhi * dx;
++iRight;
phiRight[iRight] = phi;
dPhiRight[iRight] = dPhi;
if (phiRight[iRight] * phiRight[iRight - 1] < 0)
oddNodes = !oddNodes;
}
if (phi > maxPhi)
maxPhi = phi;
if (phi < minPhi)
minPhi = phi;
}
void matchPhi () {
// divide phiLeft by maxPhi
double norm = maxPhi;
if (maxPhi < 0)
norm = -maxPhi;
maxPhi = minPhi = 0;
for (int i = 0; i <= nLeft; i++) {
phiLeft[i] /= norm;
dPhiLeft[i] /= norm;
if (phiLeft[i] > maxPhi)
maxPhi = phiLeft[i];
if (phiLeft[i] < minPhi)
minPhi = phiLeft[i];
}
// now match phiRight to phiLeft
norm = phiRight[nRight] / phiLeft[nLeft];
for (int i = 0; i <= nRight; i++) {
phiRight[i] /= norm;
dPhiRight[i] /= norm;
if (phiRight[i] > maxPhi)
maxPhi = phiRight[i];
if (phiRight[i] < minPhi)
minPhi = phiRight[i];
}
previousMisMatch = misMatch;
misMatch = dPhiRight[nRight] - dPhiLeft[nLeft];
}
void nextE () {
if (previousMisMatch * misMatch < 0) {
E -= dE;
dE /= 2;
} else {
E += dE;
}
if (Math.abs(dE) < accuracy)
eigenFound = true;
}
double Vmin = 0; // minimum energy for plotting
double Vmax = 2; // maximum energy for plotting
double xmin = -3; // minimum x for plotting
double xmax = 3; // maximum x for plotting
class Graph extends Plot {
Graph () {
setSize(400, 250);
}
public void paint () {
clear();
setColor("black");
drawAxes(xmin, xmax, Vmin, Vmax);
int points = 100;
double delta = (xmax - xmin) / points;
double xOld = xmin;
double yOld = V(xmin);
setColor("green");
for (int i = 0; i < points; i++) {
double x = xOld + delta;
double y = V(x);
plotLine(xOld, yOld, x, y);
xOld = x;
yOld = y;
}
setColor("red");
plotLine(xLeft, E, xRight, E);
}
}
class PhiGraph extends Plot {
PhiGraph () {
setSize(400, 200);
}
public void paint () {
clear();
if (maxPhi == 0 && minPhi == 0)
return;
setColor("black");
drawAxes(xmin, xmax, minPhi, maxPhi);
setColor("blue");
double dx = (xMatch - xLeft) / nLeft;
for (int i = 0; i <= iLeft; i++) {
double x = xLeft + i * dx;
plotPoint(x, phiLeft[i]);
}
dx = (xMatch - xRight) / nRight;
for (int i = 0; i <= iRight; i++) {
double x = xRight + i * dx;
plotPoint(x, phiRight[i]);
}
}
}
class Output extends Plot {
Output () {
setSize(400, 30);
}
public void paint () {
clear();
setWindow(0, 100, 0, 3);
setColor("blue");
boxArea(0, 100, 0, 3);
setColor("white");
if (iLeft > 0 && iLeft < nLeft)
plotString("integrating phi left ...", 5, 1);
else if (iRight > 0 && iRight < nRight)
plotString("integrating phi right ...", 5, 1);
else if (eigenFound)
plotString("eigen state found!!!", 5, 1);
plotString("E = " + Easy.format(E, 4), 50, 1);
plotString("dE = " + Easy.format(dE, 4), 70, 1);
}
}
Graph graph;
PhiGraph phiGraph;
Output output;
Reader E0Reader, dE0Reader, xLeftReader, xRightReader;
Reader dxReader, accuracyReader;
public void init () {
initializePhi();
add(graph = new Graph());
add(output = new Output());
add(phiGraph = new PhiGraph());
add(E0Reader = new Reader("E = ", E, 4));
add(dE0Reader = new Reader("dE = ", dE0, 4));
add(xLeftReader = new Reader("x_left = ", xLeft, 4));
add(xRightReader = new Reader("x_right = ", xRight, 4));
add(dxReader = new Reader("dx = ", dx, 4));
add(accuracyReader = new Reader("accuracy", accuracy, 4));
addControlPanel();
}
public void step () {
if (iRight < nRight) {
integrationStep();
phiGraph.repaint();
output.repaint();
} else if (!phiMatched) {
matchPhi();
phiGraph.repaint();
output.repaint();
phiMatched = true;
} else if (!eigenFound) {
nextE();
initializePhi();
graph.repaint();
output.repaint();
} else {
stopAnimation();
}
}
public void reset () {
E = E0 = E0Reader.readDouble();
dE = dE0 = dE0Reader.readDouble();
xLeft = xLeftReader.readDouble();
xRight = xRightReader.readDouble();
dx = dxReader.readDouble();
accuracy = accuracyReader.readDouble();
previousMisMatch = misMatch = 0;
phiMatched = eigenFound = false;
initializePhi();
graph.repaint();
}
public static void main (String[] args) {
Eigen eigen = new Eigen();
eigen.frame("Eigenvalue and Eigenfunction", 430, 680);
}
}
Questions or comments:
phygons@acsu.buffalo.edu