4th order runge-kutta not improving w/ step size

1. Sep 24, 2011

sam_bell

1. The problem statement, all variables and given/known data

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.

3. 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):

// 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: Sep 24, 2011