I'm modelling rocket motion using C++ and I'm solving the following differential equation:
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:
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.