Runge-Kutta method to solve Lane-Emden equation

Question

I am trying to obtain a code that will give me the numerical solutions of the Lane-Edmen equation, using the Runge-Kutta 4 steps method. I am not super into C++ and therefore there must be something I am missing, since the code is running, but the data calculated are clearly wrong. Untill now I am considering only the case n=0, for simplicity. Could you check it and give me an opinion, please?

#include <iostream>
#include <cmath>
#include <iomanip>
#include <fstream>
#define Nmax 400

using namespace std;

void RK4Step(double x, double *Y, double (*dYdx)(double, double), double (*dZdx)(double, double, double), double dx, int neq);

double dYdx(double x, double Y);

double dZdx(double x, double Y, double Z);



int main(){

    cout << setprecision(6);

    double Y[] = {1., 0.};
    double xi0 = 0.1;
    double xif = 20.;
    double dxi = (xif-xi0)/Nmax;
    int neq =  1;

    ofstream fdata;
    fdata.open("lane_emden.dat");

    fdata << xi0 << " " << Y[0] << endl;
    cout << "xi: " << scientific << xi0 << "; theta: " << Y[0] << endl;

    for(int i=0; i<Nmax; i++){
    
        RK4Step(xi0, Y, dYdx, dZdx, dxi, neq);
        xi0 += dxi;

        fdata << xi0 << " " << Y[0] << endl;
        cout << "xi: " << scientific << xi0 << "; theta: " << Y[0] << "; f: " << Y[1] << endl;

    }
    fdata.close();



    return 0;
}


void RK4Step(double xi, double *Y, double (*dYdx)(double, double), double (*dZdx)(double, double, double), double dxi, int neq){

    double k1[2], k2[2], k3[2], k4[2];
    
    k1[0] = dYdx(xi, Y[0]);
    k1[1] = dZdx(xi, Y[0], Y[1]);

    k2[0] = dYdx(xi+dxi*0.5, Y[1] + k1[1]*0.5);
    k2[1] = dZdx(xi+dxi*0.5, Y[0] + k1[0]*0.5, Y[1] + k1[1]*0.5);

    k3[0] = dYdx(xi+dxi*0.5, Y[1] + k2[1]*0.5);
    k3[1] = dZdx(xi+dxi*0.5, Y[0] + k2[0]*0.5, Y[1] + k2[1]*0.5);

    k4[0] = dYdx(xi+dxi, Y[1] + k3[1]);
    k4[1] = dZdx(xi+dxi, Y[0] + k3[0], Y[1] + k3[1]);

    for(int i=0; i<2; i++){Y[i] += (dxi/6.)*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);}
}


double dYdx(double x, double Z){

    double dY = Z;
    return dY;
}

double dZdx(double xi, double Y, double Z){

    double dZ = - 1. - (2.*Z)/xi;
    return dZ;
}

Answer 1

You will find RK4 easier to use if it is set up to use array arguments, rather than multiple separate derivative functions. Here, I've used std::valarray, but I'm well aware that many programmers don't like that.

I've set n=1. The result is more interesting than n=0. For a comparison graph, see the Wikipedia article:

https://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation

#include <iostream>
#include <fstream>
#include <iomanip>
#include <valarray>
using namespace std;

using vec = valarray<double>;

int n = 1;         // polytropic index

//======================================================================

double ipow( double x, int n ){ return n ? x * ipow( x, n - 1 ) : 1; }

//======================================================================

vec RK4step( double x, const vec &Y, vec f( double, const vec & ), double dx )
{
   vec dY1, dY2, dY3, dY4;
   dY1 = dx * f( x           , Y             );
   dY2 = dx * f( x + 0.5 * dx, Y + 0.5 * dY1 );
   dY3 = dx * f( x + 0.5 * dx, Y + 0.5 * dY2 );
   dY4 = dx * f( x +       dx, Y +       dY3 );
   return Y + ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;
}

//======================================================================

int main()
{
   vec Y = { 1.0, 0.0 };
   int Nmax = 400;
   double xi0 = 0.01, xif = 20.0;
   double dxi = ( xif - xi0 ) / Nmax;
   auto deriv = []( double xi, const vec &Y ){ return vec{ Y[1], -2.0 / xi * Y[1] - ipow( Y[0], n ) }; };

   #define FMT scientific << setw( 12 ) << setprecision( 4 )
   ofstream fdata("lane_emden.dat");
// ostream &fdata = cout;

   double xi = xi0;
   fdata << FMT << "xi" << "  " << "theta" << '\n';
   fdata << FMT <<  xi  << "  " << Y[0]    << '\n';

   for ( int i = 1; i <= Nmax; i++ )
   {
      Y = RK4step( xi, Y, deriv, dxi );
      xi += dxi;
      fdata << FMT << xi << "  " << Y[0] << '\n';
   }
}