1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Sep 24, 2011 #1
    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):

    // 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: Sep 24, 2011
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?



Similar Discussions: 4th order runge-kutta not improving w/ step size
  1. Runge kutta (Replies: 0)

Loading...