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

[sam_bell]

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

 

[KataKoniK]

Dear Sam,

Thank you for sharing your code and results. It seems like you are using a 4-point predictor-corrector scheme to solve the given differential equation. However, I noticed that you are using an RK4 scheme for the first four points. This could be the reason for the similar error values you are getting even with different numbers of grid points.

I would suggest using the same predictor-corrector scheme for all points, including the first four, to see if there is any change in the error values. Additionally, it would be helpful to check the accuracy of your code by comparing it with a known analytical solution, if available.

Another thing to consider is the accuracy of your initial conditions. The given boundary conditions are P(xi) = 1 and dP/dx(xi) = l. Make sure that these are accurately represented in your code.

I hope this helps. Let me know if you have any further questions.