Solve RK4 in C++ for Sekhar Mass White Dwarf Star

[memarf1]

Runge Kutta Method in C++ Please Help Me

This is for an RK4 second order equation.

I have run this program several times and am getting the incorrect limit for the Sekhar Mass of the White dwarf star. I was wondering if anyone has a suggestion. Ill post my Source here as well. The program seems to work and give correct values for my test equations, but its not calculating the correct limit otherwise. There is a definate limit but its about 2.15 times to high. My limit that I speak of is the mass, and it is supposed to be .54 but this program is limiting itself at 1.165. Now, I know the sekhar mass is 1.44 but the numbers .54 and 1.165 are unitless, multiplying them by 2.65 should give the sekhar mass with units. Please post suggestions.

// Runge-Kutta.cpp
//--------------------------------------------------
//A Runge-Kutta Method for solving Differential Equations
//of the form y'=f(r;x,m) ; m(r0)=0
//--------------------------------------------------

#include <iostream>
#include <iomanip>
#include <fstream>
#include <math.h>
using namespace std;

double N, X, X1, Y2, Y, H, M, R, Rf;
int (G);



double f(double r, double x, double m);
double g(double r, double x, double m);

double h(double r1, double x1, double m1);
double j(double r1, double x1, double m1);

double rungeX1(double r, double x, double m);
double rungeY1(double r1, double x1, double m1);


//Main Function
int main(double , double , double)
{
double x, m, r, x1, m1, r1;

ofstream outfile;
outfile.open("outfile.doc");

while (G != 2)
{

cout << "\n";
cout << "\n";
cout << "\n";
//cout << "Please enter a value for R0: ";
//cin >> R;
cout << "\n";
//cout << "Please enter a value for M0: ";
//cin >> M;
cout << "\n";
cout << "Please enter a value for H: ";
cin >> H;
cout << "\n";
//cout << "Please enter a value for Rfinal: ";
//cin >> Rf;
cout << "\n";
cout << "Please enter a value for X0: ";
cin >> X;
cout << "\n";
cout << "\n";
cout << "\n";

outfile<<"\t******************* Runge-Kutta 4 *******************"
<<"\n\n";
outfile<<" "
<<setw(17)<<"r"<<setw(17)<<"x"<<setw(16)<<"\tm"
<<"\n"
<<"\t-----------------------------------------------------"
<<"\n";

cout<<"\t******************* Runge-Kutta 4 *******************"
<<"\n\n";
cout<<" "
<<setw(17)<<"r"<<setw(17)<<"x"<<setw(16)<<"\tm"
<<"\n"
<<"\t-----------------------------------------------------"
<<"\n";

//**********************************************************
M = 0;
R = 1 * pow(10, -23);
Rf = 20;
//**********************************************************

m = M;
m1 = M;
x = X;
x1 = X;
N = Rf / H;
r=R;



for(int i=0;i<=N;i++)
if(x>0)
{
r=R+(i*H) + H;
r1=R+(i*H) + H;
x=rungeX1(r,x,m);
x1=x;
m=rungeY1(r1,x1,m1);
m1=m;



outfile<<left<<setw(6)<<i<<"|"
<<setprecision(12)<<left<<setw(8)<<"\t"<<r
<<setprecision(12)<<left<<setw(8)<<"\t"<<x
<<setprecision(12)<<left<<setw(8)<<"\t"<<m;
outfile<<"\n";

cout<<left<<setw(6)<<i<<"|"
<<setprecision(12)<<left<<setw(8)<<"\t"<<r
<<setprecision(12)<<left<<setw(8)<<"\t"<<x
<<setprecision(12)<<left<<setw(8)<<"\t"<<m;
cout<<"\n\n";


}

cout<< endl<<endl<<endl<<endl;
cout<< " **************** EVALUATE YOUR DATA ****************";
cout<< endl<<endl<<endl<<endl;
cin >> G;
cout<<"\n\n";

outfile.close();
return 0;
}


}


double rungeX1(double r, double x, double m)
{
double K1 = (H * f(r,x,m));
double M1 = (H * g(r,x,m));

double K2 = (H * f((r + H/2),(x + (K1)/2),(m + (M1)/2)));
double M2 = (H * g((r + H/2),(x + (K1)/2),(m + (M1)/2)));

double K3 = (H * f((r + H/2),(x + (K2)/2),(m + (M2)/2)));
double M3 = (H * g((r + H/2),(x + (K2)/2),(m + (M2)/2)));

double K4 = (H * f((r + H),(x + (K3)),(m + (M3))));
double M4 = (H * g((r + H),(x + (K3)),(m + (M3))));

double X1 = ((((K1) + (2 * K2) + (2 * K3) + (K4)))/6 + x);
double Y1 = ((((M1) + (2 * M2) + (2 * M3) + (M4)))/6 + m);

//cout << K1 << ' ' << K2 << ' ' << K3 << ' ' << K4 << ' ' << X1 << endl;
//cout << M1 << ' ' << M2 << ' ' << M3 << ' ' << M4 << ' ' << endl;

return X1;
}

double rungeY1(double r1, double x1, double m1)
{
double P1 = (H * h(r1,x1,m1));
double L1 = (H * j(r1,x1,m1));

double P2 = (H * h((r1 + H/2),(x1 + (P1)/2),(m1 + (L1)/2)));
double L2 = (H * j((r1 + H/2),(x1 + (P1)/2),(m1 + (L1)/2)));

double P3 = (H * h((r1 + H/2),(x1 + (P2)/2),(m1 + (L2)/2)));
double L3 = (H * j((r1 + H/2),(x1 + (P2)/2),(m1 + (L2)/2)));

double P4 = (H * h((r1 + H),(x1 + (P3)),(m1 + (L3))));
double L4 = (H * j((r1 + H),(x1 + (P3)),(m1 + (L3))));

double X2 = ((((P1) + (2 * P2) + (2 * P3) + (P4)))/6 + x1);
double Y2 = ((((L1) + (2 * L2) + (2 * L3) + (L4)))/6 + m1);

//cout << P1 << ' ' << P2 << ' ' << P3 << ' ' << P4 << ' ' << endl;
//cout << L1 << ' ' << L2 << ' ' << L3 << ' ' << L4 << ' ' << Y2 << endl;

return Y2;
}





double f(double r, double x, double m)
{
double f = ((-5/3) * (1 / x) * (m / pow(r, 2)) * (sqrt(1 + pow(x, 2))));
//double f = -r + x + m;
//double f = 2*x + 4*m;
return f;
}

double g(double r, double x, double m)
{
double g = 3 * pow(r, 2) * pow(x, 3);
//double g = 2 * r + 3 * x - m;
//double g = -x +6*m;
return g;
}

double h(double r1, double x1, double m1)
{
double h = ((-5/3) * (1 / x1) * (m1 / pow(r1, 2)) * (sqrt(1 + pow(x1, 2))));
//double h = -r1 + x1 + m1;
//double h = 2*x1 + 4*m1;
return h;
}

double j(double r1, double x1, double m1)
{
double j = 3 * pow(r1, 2) * pow(x1, 3);
//double j = 2 * r1 + 3 * x1 - m1;
//double j = -x1 +6*m1;
return j;
}

 

[memarf1]

Please help me, does anyone know what is wrong with this?

 

[M98Ranger]

Hello,

Thank you for reaching out for assistance with your RK4 program in C++. I have reviewed your code and I have a few suggestions that may help you solve the issue you are experiencing.

1. Double check your initial conditions: In your main function, you have set the initial values of R and M to be very small (1e-23 and 0, respectively). These values may not be appropriate for your problem and could be causing your program to give incorrect results. I suggest double checking your initial conditions and making sure they are appropriate for the problem you are trying to solve.

2. Check your equations: In your functions f and g, you have used the values -5/3 and 1/2, respectively. These values may not be the correct coefficients for your problem and could be contributing to your incorrect results. I suggest double checking your equations and making sure they are correct.

3. Use more precise data types: In your main function, you have declared all your variables as type double. While this may be appropriate for some problems, it may not be precise enough for your problem. I suggest using a more precise data type, such as long double or even a library that supports arbitrary precision, like Boost.

4. Use a smaller step size: In your main function, you have set the step size H to be a user input. Depending on the problem, this step size may be too large and could be causing your program to give incorrect results. I suggest trying a smaller step size and seeing if that helps improve the accuracy of your results.

I hope these suggestions help you solve the issue you are experiencing with your RK4 program. If you continue to have trouble, I suggest seeking help from a tutor or mentor who has experience with solving RK4 equations in C++. Good luck!