Hi there
I'm trying to solve the rocket motion equation in C using the Numerical Recipes engine. I just can't get the right output-results - they're negative.
I assume that the exhaust velocity and rate of mass change is constant.

This is the equation:

http://it-drengene.dk/1.jpg

where ve is the exhaust velocity, alpha is the burn rate, M0 is the mass of rocket + fuel, g is the gravitational constant, c is the drag constant and vr is the velocity of
the rocket.

I think my error is in the calculations.

This i what i've got:


/*
Equation file
*/
#include <math.h>
#include "funcs.h" /* contains number of parameters ... 6 in this case */

#define ALPHA eqpar[0] /*burn rate */
#define M_0 eqpar[1] /* mass of rocket + fuel */
#define g eqpar[2] /* gravitational constant */
#define C eqpar[3] /* drag constant */
#define V_e eqpar[4] /* exhaust velocity */
#define M_t eqpar[5] /* mass of empty rocket */

void eq(double t, double y[], double dydt[])
{
extern double eqpar[NPAR];
double gravity, drag, rocket, deltamass;

deltamass = M_0 - (ALPHA * t);
if(deltamass<M_t) { deltamass = M_t; }

drag = (C * y[2]*y[2])/deltamass;
if(y[2] > 0) { drag= -1; }

gravity = -g;

rocket = (V_e*ALPHA)/deltamass;
if (deltamass = M_t) { rocket =0; }

dydt[1] = y[2];
dydt[2] = gravity + drag + rocket;
return;
}

#undef ALPHA
#undef M_0
#undef g
#undef C
#undef V_e
#undef M_t


and the dat-file containing parameters:


datafile containing parameters

# 1 3 data points
# 1 ALPHA
# 10 M_0
# 9.82 g
# 0.001735 C
# -5 V_e
# 5 M_t
0.000000e00 1.000000e00 -1.00000000


This is my output when taking steps of 0.1 seconds with 10 steps:


# 10 3 data points
# 1.000000e+000
# 1.000000e+001

time - pos - velocity
1.0000000000e-002 -2.4382614634e-020 -2.4382614634e-018
2.0000000000e-002 -4.8765229268e-020 -2.4382614634e-018
3.0000000000e-002 -7.3147843902e-020 -2.4382614634e-018
4.0000000000e-002 -9.7530458537e-020 -2.4382614634e-018
5.0000000000e-002 -1.2191307317e-019 -2.4382614634e-018
6.0000000000e-002 -1.4629568780e-019 -2.4382614634e-018
7.0000000000e-002 -1.7067830244e-019 -2.4382614634e-018
8.0000000000e-002 -1.9506091707e-019 -2.4382614634e-018
9.0000000000e-002 -2.1944353171e-019 -2.4382614634e-018
1.0000000000e-001 -2.4382614634e-019 -2.4382614634e-018


As you can see, my the output is different from what you could expect. As i said, i think the error lies in the equation file.

This is my output when taking steps of 0.1 seconds with 10 steps:

Your data looks like your steps are in increments of 0.01 not 0.1 (1.0e-2 = 0.01).

From your output it looks like when the time = 0.01 seconds the results for the position and velocity are zero, as expected.

Why don't you integrate furthur, say to 20 seconds, and plot the results to see what the position and velocity are doing. Unless this rocket takes off and within 0.01 seconds accelerates to near the velocity of light, not much is going to happen in .1 seconds.

I haven't studied your equations in detail but ....

Your gravitational constant is a positive 9.82 m/s^2. Have you tried changing it to a negative 9.82 m/s^2?

.... output is different from what you could expect

What do you expect? Do you have a hand calculation for comparison of results?

Thanks
Matt

Your data looks like your steps are in increments of 0.01 not 0.1 (1.0e-2 = 0.01).


Yes, your right - i used 0.01 :) My mistake.


Why don't you integrate furthur, say to 20 seconds, and plot the results to see what the position and velocity are doing. Unless this rocket takes off and within 0.01 seconds accelerates to near the velocity of light, not much is going to happen in .1 seconds.


Yes, i have plotted the graph
http://it-drengene.dk/2.jpg



Your gravitational constant is a positive 9.82 m/s^2. Have you tried changing it to a negative 9.82 m/s^2?



Look at the sourcecode
gravity = -g;

I can expect the velocity as a function of time - and that should be positive :)