- #1
memarf1
- 18
- 0
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;
}
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;
}
Last edited: