**where art meets electronics**

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” 2^{nd}
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 ρ:

Where

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 4^{th}-order Runge-Kutta integration. I recommend you
watch this video to
learn about it.

For a 1D ODE the
Runge-Kutta 4^{th}-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·10^{8} and 1·10^{16}
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 (christoph@christophstark.de)

#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;

result.dz = -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*k1.dz*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*k2.dz*dStepSize;

Derivative k3= DGL_RHS(val);

val.x = state.x + dStepSize;

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

val.z = state.z + k3.dz*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);

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

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

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

{

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

return;

}

//Advance the state

state.x += dStepSize;

state.y += deriv.dy*dStepSize;

state.z += deriv.dz*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

dataFile.open(s,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

{

i=0;

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;

//OutputDebugString(s);

}

i++;

//Do one integration step

RK4_integration(state,dr);

}

dataFile.close();

//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

printf(s);

getchar(); //wait for user press enter

}

Last change: 2010-11-08