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.

 

[DeeNos]

I am impressed by your use of the Runge-Kutta method to study the structure of white dwarf stars. However, I noticed that you have used a non-relativistic approximation in your code. While this may be a good starting point, I would suggest considering the effects of relativity in your calculations as they may play a significant role in the structure of white dwarf stars. Additionally, have you considered incorporating other important factors such as temperature and composition into your equations? These can greatly affect the behavior of white dwarf stars and should be taken into account in your calculations. I would also recommend running your code repeatedly to see if the results converge to more reasonable values. Overall, your approach seems promising and I look forward to seeing the results of your calculations. Best of luck in your research!