- #1
sam_bell
- 67
- 0
Homework Statement
Hi. I've written a program to find the solution of the differential equation
[ -d/dx^2 + d/dx + l(l+1) + r^2(V(x)-E) ]P(x) = 0
where l is an integer constant, V(x) is defined on a grid, and E is a real number.
The boundary condition is P(xi) = 1, dP/dx(xi) = l.
I'm using a 4 point predictor-corrector scheme for most of the integration except that it needs four points to get started. For that I'm using an RK4 scheme from the initial position. To estimate the error of my solution, I printed the LHS of the above equation. (To get the derivative a local, cubic fit is used.) What I don't understand is that this error for the first four points is the same even as I go from 1000 to 10000 to 100000 grid points.
The Attempt at a Solution
I changed the above equation into a 1st order equation:
d/dx [P(x),Q(x)] = [Q(x), ( l(l+1) + r^2(V(x)-E) )P(x) + Q(x)]
Here is how I implemented it (i think it's readable):
// perform 4-step Adams-Bashforth-Moulton integration
// change independent variable to x = log(r/mu):
// [-(d/dx)^2 + d/dx + l(l+1) + r^2 (V-E)] P = 0
// define Q = dP/dx to change to first order equation:
// d[P, Q]/dx = [Q, fP+Q]
// where f(x) = l(l+1) + r(x)^2 ( V(x)-E )
double fP[4], al[4], c[4], dP, dQ;
double k1P,k1Q,k2P,k2Q,k3P,k3Q,k4P,k4Q,f23;
// precompute f(x) for all x
for(j = 0; j < N; j++) f[j] = C + r[j]*(rV[j]-r[j]*E);
// get cubic fit of f(x) for first four grid points
al[0] = xi; al[1] = xi+dx; al[2] = xi+2.0*dx; al[3] = xi+3.0*dx;
getcubicfit(al, f, c);
// set boundary condition P_nl(r=0) = 0, dP/dr(r=0) = A
// recall index=0 corresponds to r=r_min, and Q = dP/dx
P[0] = 1.0; Q[0] = l*1.0; fP[3] = f[0]*P[0];
// compute first four values of [P,Q] using RK4
x = xi;
for(j = 1; j < 4; j++) {
// use cubic fit to get f(x+dx/2)
for(f23 = 0.0, k = 3; k >= 0; k--)
f23 += f23*(x+dx/2.0) + c[k];
k1P = dx*Q[j-1];
k1Q = dx*f[j-1]*P[j-1] + k1P;
k2P = dx*(Q[j-1]+k1Q/2.0);
k2Q = dx*f23*(P[j-1]+k1P/2.0) + k2P;
k3P = dx*(Q[j-1]+k2Q/2.0);
k3Q = dx*f23*(P[j-1]+k2P/2.0) + k3P;
k4P = dx*(Q[j-1]+k3Q);
k4Q = dx*f[j]*(P[j-1]+k3P) + k4P;
P[j] = P[j-1] + (k1P + 2.0*k2P + 2.0*k3P + k4P)/6.0;
Q[j] = Q[j-1] + (k1Q + 2.0*k2Q + 2.0*k3Q + k4Q)/6.0;
fP[3-j] = f[j]*P[j]; x + dx;
}
And here is the output for the first few points
100000 grid points
Indx x P(x) LHS
2 -6.907511155087987 0.999999997976578 0.035358035254560
3 -6.907389093140912 0.999999995447077 0.022094527146322
4 -6.907267031193837 0.999999992103051 -0.005888350745711
5 -6.907144969246762 0.999999988361392 0.001473385819706
6 -6.907022907299686 0.999999984112253 -0.000002117234924
7 -6.906900845352611 0.999999979377458 -0.000001373001447
8 -6.906778783405536 0.999999974156891 0.000000646808609
10000 grid points
Indx x P(x) LHS
2 -6.905313820307163 0.999999797345727 0.035348203876479
3 -6.904093090969677 0.999999543803950 0.022084443517820
4 -6.902872361632190 0.999999208328583 -0.005890165535733
5 -6.901651632294704 0.999998832468881 0.001471261077949
6 -6.900430902957218 0.999998405150760 0.000000309298519
7 -6.899210173619731 0.999997928443006 -0.000000063535007
8 -6.897989444282245 0.999997402227227 0.000000009280229 1000 grid points
Indx x P(x) LHS
2 -6.883318697109203 0.999979432800277 0.035211195309530
3 -6.871100406172737 0.999953495038619 0.021923331632628
4 -6.858882115236270 0.999918874396041 -0.005853988700196
5 -6.846663824299803 0.999879570436276 0.001452321145710
6 -6.834445533363336 0.999834376636500 0.000004247349503
7 -6.822227242426870 0.999783372570865 0.000000589394914
8 -6.810008951490403 0.999726424184519 0.000001334948859
Getting real confused. Thanks for any suggestions.
Sam
Last edited: