Numerical Integration of two coupled ordinary first-order differential equations

Figure 1: Calculated mass and density inside a carbon white dwarf as function of position.

Sometimes you come across a problem that has no obvious analytic solution. In these cases a numerical solution can be the answer. Here I want to discuss an example problem (6.11) from the book “The Physics of Stars” 2nd Edition by A.C. Philips, where we have to find the internal structure of white dwarf stars.

The mass and density distribution in a star is governed by the following two coupled ordinary first-order differential equations (ODEs):

These equations need to be solved simultaneously, which makes it a little tricky. We are also given the derivative of the pressure with regard to ρ:


is the dimensionless Fermi moment, is the hydrogen mass, and is the average number of electrons per nucleon.

The Idea is to calculate the right hand side of the ODE and then use this value to calculate the change due to a small step. There are different ways of doing this, the easiest being an Euler integration. This method can easily give inaccurate results and therefore the common method nowadays is the 4th-order Runge-Kutta integration. I recommend you watch this video to learn about it.

For a 1D ODE   the Runge-Kutta 4th-order (RK4) formulas are:

where h is the step size,  and  are the current values, and  and  are the values at the future step. For more than 1 dimension we treat  as a vector.

To implement the formulas start by writing our equations in dimensionless form. We write the radius in meters, mass in sun masses, and density in units of the central density.

Our ODEs look like this now:

We use the RK4 formula to advance the state of our variables by small steps. We start with the boundary condition that the central density at x=0 is given. We the work our way outwards until the density falls below a certain threshold. The tricky part is to specify the behavior at x=0. If we really start at x=0, the first ODE has a divide by zero, and the second ODE is zero. The easy solution is to start a little away from x=0. White dwarfs have radii of thousands of km, so it won’t hurt if we start a few micrometer away from the origin, e.g. at x=1·10-5. We also have to set an initial mass, so that the first iterations of the calculation don’t blow up. The following code is tested with MS VS 2008 (compile as console application) and saves the mass and density distribution to a text file so you can easily plot them (for example with the free QtiPlot). The values of the radius and the mass are printed to the screen. I hope I have added enough comments to the code so that it’s easy to understand.

The function DGL_RHS calculates the expression at the right hand side of our dimensionless ODEs and returns the derivatives in a structure. It needs the values of the current state (for x,z, and y) as input.

The function RK4_integration(Values &state, double dStepSize) calculates the new values after one single step. The new values will be saved back to the variable called “state”.

The main() function runs a loop that calls RK4_integration. It stops when the density falls below a threshold. This means that we reached the surface of the star, because less and less material is above of the current position.

We can vary the central density from 1·108 and 1·1016 kgm-3 and plot the resulting masses vs. the radius.

Figure 2: The radius and mass of white dwarfs are correlated. Here are some simulation data for carbon and iron dwarfs.

// WhiteDwarf RK4 integration.cpp

// Copyright (c) 2010, Christoph Stark (


#include <math.h>

#include <float.h>

#include <Windows.h>

#include <iostream>

#include <fstream>

using namespace std;


//all constants in this program are in SI units!

const double PI=3.1415926;

const double G=6.673E-11;

const double Y_e_iron = 26.0/56.0; //26/56 for Iron white dwarf

const double Y_e_carbon = 6.0/12.0; //6/12 for carbon white dwarf

const double Y_e = Y_e_carbon;

const double m_e=9.109E-31; //kg

const double m_H=1.673E-27; //kg, Proton

const double c = 299792458; //m/s

const double h = 6.6261E-34;

const double M_sun =1.989E30; //kg

double rho_c = 0; //kgm^-3 will be varied


struct Values


      double x; // Radius in meters

      double y; // mass [M_sun]

      double z; // density [rho_c]



struct Derivative


      double dy;

      double dz;



//calculate the right side of the DGL = derivatives

Derivative DGL_RHS(Values &start)


      Derivative result;

      //calc. dimension less Fermi momentum xF

      double x_F=pow(3*Y_e*start.z*rho_c/(8.0*PI*m_H),1.0/3.0)*h/(m_e*c);

      //calculate dRho/dP fist

      double dP_drho = Y_e*(m_e*c*c/m_H)/3.0*(x_F*x_F/sqrt(1.0+x_F*x_F));

      result.dy = 4.0*PI/M_sun*rho_c*start.x*start.x*start.z; = -G*M_sun*start.y*start.z/start.x/start.x/dP_drho;

      return result;



void RK4_integration(Values &state, double dStepSize)


      //Calculate the RK4 coefficients

      Values val;

      Derivative k1= DGL_RHS(state);

      val.x = state.x + 0.5*dStepSize;

      val.y = state.y + 0.5*k1.dy*dStepSize;

      val.z = state.z + 0.5**dStepSize;

      Derivative k2= DGL_RHS(val);

      val.x = state.x + 0.5*dStepSize;

      val.y = state.y + 0.5*k2.dy*dStepSize;

      val.z = state.z + 0.5**dStepSize;

      Derivative k3= DGL_RHS(val);

      val.x = state.x + dStepSize;

      val.y = state.y + k3.dy*dStepSize;

      val.z = state.z +*dStepSize;

      Derivative k4= DGL_RHS(val);


      //RK4 weights

      Derivative deriv;

      deriv.dy = 1.0/6.0*(k1.dy+2.0*(k2.dy + k3.dy) + k4.dy); = 1.0/6.0*(*( + +;

      //if we get invalid numbers, we probably reached the surface of the star

      if ((_isnan(deriv.dy ))||(_isnan(


            state.z=0; //signal to quit the main loop. state.y will contain last valid calculation



      //Advance the state

      state.x += dStepSize;

      state.y += deriv.dy*dStepSize;

      state.z +=*dStepSize;



int main(int argc, char* argv[])


      rho_c = 1E10; // chose a central density


     //start with a new star

      Values state;

      state.x = 1E-5; //avoid r=0 at center of star. This is somewhat critical. If too small, the first derivative will blow up!

      state.y = 1E-5/M_sun; //start with a non-zero mass. To avoid nummerical problems at beginning

      state.z = 1; //1=one times rho_c

      double dr = 1;


      int len = 100;

      char* s=new char[len];

      int i=0;


      //save data to a file, so we can plot it later

      std::ofstream dataFile;

      sprintf_s(s,sizeof(char)*len,"White dwarf rho_c=%G.csv",rho_c); //make a filename,ios_base::out | ios_base::trunc); //overwrite existing file

      dataFile << "Radius\tMass\tRho\n";

      dataFile << "m\t[M_sun]\t[rho_c]\n";


      //run the simulation until the density drops below 1E-9 kgm^-3 at surface of star

      const double stop_z = 1E-9;

      while (fabs(state.z)> stop_z )


            if (5000 == i) //we only save a few data points



                  sprintf_s(s,sizeof(char)*len,"%10.5E \t %10.5E\t %10.5E\n",state.x, state.y, state.z);

                  //save data to file

                  dataFile << s;




            //Do one integration step





      //print final value of this white dwarf

      sprintf_s(s,sizeof(char)*len,"%10.06E \t %10.6E \t %10.4E\t %10.4E\n",rho_c,state.x, state.y, state.z);

      //print the values for this central density

      //be aware of dimensionless units



      getchar(); //wait for user press enter




Last change: 2010-11-08