Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Runge Kutta Method

  1. Aug 5, 2005 #1
    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;
    }
     
    Last edited: Aug 5, 2005
  2. jcsd
  3. Aug 7, 2005 #2
    Please help me, does anyone know what is wrong with this?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Runge Kutta Method
  1. Runge Kutta Method (Replies: 3)

Loading...