White Dwarf Star Coupled Runge-Kutta

[snelson989]

Homework Statement

Use a coupled fourth order Runge-Kutta, to find the structure of white dwarf stars.
I think I am applying the Runge-kutta method wrong?
Variables defined in C code notes

Homework Equations

Equations in c code. and in attached images.

The Attempt at a Solution

This is the code for running through the 4th order runge kutte once, I have not done this repeatedly yet, as the numbers I get for the first run are ridiculous.
I have used a non relativistic approximation.

double rho0, rho1, rho2, rho3, rho4, f0, f1, f2, f3, h=6.626068*pow(10,-34), me=9.10938188*pow(10,-31), mp=1.67262158*pow(10,-27), pi=3.14159265,
r0, r1, r2, r3, g0, g1, g2, g3, dr=1, rhoc, G=6.67300*pow(10,-11), m0, m1, m2, m3, m4;
/* rho# values of density for Runge-Kutta method.
f# values of the derivative of pressure with respect to density used in Runge-Kutta method.
me electron mass
mp photon mass
h Plancks constant
r# radii for Runge-Kutta method
g# derivative of mass with respect to radius for Runge-Kutta method.
dr increase in radius with Runge-Kutta step
drho increase in density with Runge-Kutta step.
rhoc central density
m# mass contained in the star integrated up to that point
G gravitational constant
*/

for( r0=pow(10,-10), rho0=pow(10,10), m0=(4/3)*pi*pow(10,-20); rho0>=0; rho0=rho4, r0=r3)
{
g0=4*pi*r0*r0*rho0;
f0=(-1*(48*me*(pow(mp,(5/3))))/(h*h*(pow(2,1/3))*(pow((3*rho0/pi),(2/3)))))*G*rho0*m0/(r0*r0);

r1=r0+dr/2;
r2=r0+dr/2;
r3=r0+dr;

m1=m0+g0*dr/2;
rho1=rho0+f0*dr/2;
g1=4*pi*r1*r1*rho1;
f1=(-1/(h*h*(pow(2,1/3))*(pow((3*rho1/pi),(2/3)))/(48*me*(pow(mp,(5/3))))))*G*rho1*m1/(r1*r1);

m2=m0+g1*dr/2;
rho2=rho0+f1*dr/2;
g2=4*pi*r2*r2*rho2;
f2=(-1/(h*h*(pow(2,1/3))*(pow((3*rho2/pi),(2/3)))/(48*me*(pow(mp,(5/3))))))*G*rho2*m2/(r2*r2);


m3=m0+g2*dr/2;
rho3=rho0+f2*dr/2;
g3=4*pi*r3*r3*rho2;
f3=(-1/(h*h*(pow(2,1/3))*(pow((3*rho3/pi),(2/3)))/(48*me*(pow(mp,(5/3))))))*G*rho3*m3/(r3*r3);

m4=m0+(g0+2*g1+2*g2+g3)*dr/6;
rho4=rho0+(f0+2*f1+2*f2+f3)*dr/6;
printf("%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n", rho0, rho1, rho2, rho3, rho4, f0, f1, f2, f3, r0, r1, r2, r3, g0, g1, g2, g3, m0, m1, m2, m3, m4);
system("PAUSE");
}

fprintf(fout,"%f\t%f\t%f\n", m4, rho4, r3);
system("PAUSE");
exit(0);
}

 

[D H]

First off, you are not writing in C. You are writing in C-tran.

That, however, is not your problem.

One problem is that 2/3 = 1/3 = 0.

Another is that your expression for f0 looks very different from your expressions for f1, f2, and f3.

 

[snelson989]

First off, you are not writing in C. You are writing in C-tran.

That, however, is not your problem.

One problem is that 2/3 = 1/3 = 0.

Another is that your expression for f0 looks very different from your expressions for f1, f2, and f3.

thankyou, I assume that I need to define the fractions as floats or doubles.
however I think that although the functions for f0 etc look different I have simply swapped a pow(x,-1) for explicitly doing that operation, although I may well be wrong I have changed many functions in an a ttempt to isolate the problem.

Thankyou You've finally solve a weeks anguish of me struggling with this.

 

[D H]

You're welcome.

It is your willy-nilly usage of pow() that made me say you are writing C-tran. That C has a pow() function as opposed to an exponential operator is one of the key weaknesses of C/C++ when it comes to scientific computing. pow(a,b) is typically evaluated as exp(b*ln(a)). The function is something to avoid if at all possible.