Numerican Int of ODE

[Dulle]

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.

[waterborne]

There are a few places you could be messing up (or I could be mistaken ). You are trying to solve this ODE but I don't see any numerical solver code. If you don't have a numerical solver that could be part of the problem. What I see in your code is just values of the derivative of velocity. Check out the Runge-Kutta method for solving ODEs. If you are in a time crunch Euler's method is easier (and faster to code) but not as accurate.

My second issue is with this line

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

I am not an aerospace engineer but I thought drag was generally proportional to v not v^2. Regardless of whether or not that is correct, it looks like when you have an increasing velocity your drag remains constant. Intuitively, that does not make sense. I could be misreading your code and not fully understand your notation.

Lastly, my you could try setting g=0 and drag=0 for all time. Then you know you should absolutely get positive velocity. If you don't then either there is a sign error somewhere in your code and/or a programming mistake. Hope this helps