Solving a second order ODE using a Runge-Kutta algorithm in C++

1.5k views Asked by At

I'm modelling rocket motion using C++ and I'm solving the following differential equation:

Second order ODE.

I want to solve this using Runge-Kutta and to do that I'm going to use a gsl routine. First I need to decompose it into two first order ODEs which results in:

First order ODEs

Where tau is the only variable (U = 1 in this case). Note that tau = t/tb (this is important since I enter it as t/tb in my code)! So I'm going to use a Runge-Kutta 4 algorithm, the problem is, I think I'm making a mistake when inputting these two equations into vector form into the code. If somebody could spot what i'm doing wrong I would be forever grateful.

#include <stdio.h>
#include <gsl/gsl_erno.h>
#include <gsl/gsl_matrix.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>

const double tg = 804.9;
const double tc = 10000;
const double g = 9.829;
const double U = 1;
const double r0 = 0.6368*pow(10,7);
const double tb = 200;

int func(double t, const double y[], double f[], void *params)
{
  f[0] = y[1];
  f[1] = tb*tb/(tg*tg)*(tc*U/(tb*(1-t/tb)) - 1);
  return GSL_SUCCESS:
}

int main (void);
{
  double A;

  printf( "%s\n", "Please enter the desired step size" );
  std::cin >> A;


  gsl_odeiv2_system = {func, 0, 2, &mu};


  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys,gsl_odeiv2_step_rk4, 1e-3, 1e-8, 1e-8);

  double t = 0.0;
  double y[2] = {r0, 0.0};
  printf( "%f $10f\n", y[0], y[1] );
  int i, s;

  printf( "%s %22s %22s\n", "Time", "x(t)", "v(t)" );

  for (i = 0; i <= tb; i++)
  {
    printf( "%.15f %22.15f %22.15f\n", t/tb, y[0]/r0, y[1]*tb/r0 );

    s = gsl_odeiv2_driver_apply_fixed_step (d, &t, A, 1/A, y);

  }
  gsl_odeiv2_driver_free (d);

return 0;
}

I posted my whole code just in case. Keep in mind that my code says t/tb but t/tb = tau in fact. All the constants makes it look confusing and messy but the bit that I think I have gotten wrong is in the function 'func'. So I suppose my question is; by looking at the first order ODEs I've uploaded, have I inputted them into my code correctly? Thank you very much in advance.

0

There are 0 answers