Hi there
Im currently working at a water-rocket simulation. I'm having trouble writing the right code for eq.c.
The water-rocket equation is as follows:

155

and divided by dt the rocket acceleration in accelerationinterval

156

Where dvr is the velocitychange of the rocket, ve is the velocity of the water ejected from the rocket, Mr is the mass of the rocket and dMv is the changing mass of the water.

I figured eq.c should look something like this:

#include <math.h>
#include "funcs.h"

#define ... eqpar[0]

void eq(double t, double y[], double dydt[])
{
extern double eqpar[NPAR];
dydt[1] = ...
return;
}

#undef ...

If you are looking for help, here is the methodology you should use.

From Davekw7x (with modifications to provide a generalized version of the methodology).

This quote, from Davekw7x, is taken from this thread Chap.10: Powell minimization post #5.

Here's the drill:

1. Tell us what the program is supposed to do.

2. For purposes of getting help, I suggest that you reduce the program to the part that is problematic.

3. Show us exactly what you expect to see in the output. Of course you might not know in advance what the exact output values are, but show us what you expect the output to look like.

4. Show us exactly what you got from the program.

5. Tell us what exactly you don't understand about the output.

6. Make the output on your program verbose enough so that you (and we) can know exactly what the program is telling us

We can't just help you without knowing, at a minimum, the above information.

Thanks
Matt

Sorry, i wanted to model the rocket in a single dimension, e.g. the velocity of the rocket over time dt. So i need to solve the differential equation and model it's acceleration/velocity over time dt and graph it in gnuplot.

...working at a water-rocket simulation. ...

Hmmm...

Your first attachment appears to be the basic rocket equation with no external forces: (See Footnote.)


\mathrm{d}v = -V_e \frac{\mathrm{d}M}{M}


Where:

v
is the velocity of the rocket with respect to ground.


-V_e
is the exhaust velocity of the water with respect to the rocket.


M
is the total mass (mass of the empty rocket plus the mass of the water).


\mathrm{d}v
is the change in velocity over an interval \mathrm{d}t


\mathrm{d}M
is the change in mass over the interval \mathrm{d}t



If we start with an initial velocity of v_0 and an initial mass of M_0, then, if V_e is constant, we can integrate the above equation to get

v(t) = v_0+V_e\, \ln \frac{M_0}{M(t)}

Now: Where would you like to go from here? In other words, what exactly are you given, and what exactly do you want to solve for?


Regards,

Dave


Footnote: Where did you get the equation for your second attachment, and what's the point of it? I mean, you can't just take an equation involving variables and differentials and "divide by" \mathrm{d}t as though it's an algebraic expression.

There's an interesting and fairly comprehensive treatment of the subject here: http://www.sciencebits.com/RocketEqs. It starts with fundamentals leading up to your first equation, then puts in the external forces (gravity, drag) to get to a set of equations that may give some practical results.

I simply want to model the rocket-equation in C using NR - that's why i made eq.c - this file contains the equation to solve and model.
I just need to "write" the rocket equation in eq.c to solve it. Sorry for not being understandable.

I'm not interested in adding external forces before i get this thing solved.

I simply want to model the rocket-equation in C...get this thing solved.
I hate to repeat myself, but what "thing" do you want to solve? We have an equation that shows v(t) as a function of initial velocity, initial mass, nozzle velocity and mass as a function of time, (derived with certain assumptions about external forces, etc.).

Are you wanting to integrate v(t) to find distance as a function of time or what?

In order to do anything numerically, you have to have some numbers. In order to program some function, you have to define the function in a way that can be evaluated numerically.

So:

What is the initial velocity? I'm thinking that you have it equal to zero. That's not too tough.

What is the Initial mass? Give us a value to work with.

What is the nozzle velocity? Is it equal to some constant during the thrust phase? My function for v(t) obtained from dv assumes that it is, but, in general, it doesn't have to be. If you want to obtain a function for v(t) from the equation for dv, you have to tell us something about it. What is its value? If it is not constant then define it somehow. Be specific. Give us something to work with. (Some number to plug into the program.)

What is the total mass? We know that it is a function of time, since the amount of water is decreasing throughout the thrust phase. What kind of function is it? Is it linearly decreasing during the thrust phase or what? Define it exactly so that we have something to work with. Define it in words or in some mathematical terms. Give us some numbers.

Your question is about how to use NR for something. My question is about what problem you are trying to solve and exactly what you (we) are given.


Regards,

Dave

I've got the basic rocket equation with no external forces.
I want to model my rocket in a single dimension - i want to find the velocity as a function of time. So i figures i need a default.dat file with three "columns" - time, position and velocity.

Let's say that the change in mass -dM and the exhaust velocity of the water with respect to the rocket is constant.
So i need to define the initial velocity and mass.

The default.dat-file should probably look something like this:

# 1 3 data points
# 0 VO
# 1.1 MO
0.000000e00 1.000000e00 -1.00000000

where VO is the initial velocity and MO is the initial mass of the rocket. The three columns define the output - time, position and velocity and VO and MO is the input-parameters.
Now i want to "write" the rocket-equation in a file eq.c containing the equation and equation-parameters. THIS is what i'm unsure of - to write the differential-equation in c.
I figured that eq.c should look something like this:

#include <math.h>
#include "funcs.h"

#define VO eqpar[0]
#define MO eqpar[1]

void eq(double t, double y[], double dydt[])
{
extern double eqpar[NPAR];
dydt[1] = y[2];
dydt[2] = ...; /* equation */
return;
}

#undef VO
#undef MO

The funcs.h file simply contains the number of parameters etc.
Anyway - i'm very grateful for your time.

... want to find the velocity as a function of time.If you are given a first order ordinary differential equation in the form dv/dt = f(t) and you are given v(0), you can try integrating it numerically with something like Simpson's rule or some such thing.

Well, here goes...

To solve the rocket equation for v(t), that we saw earlier, you have to be given the following:

v(0): The initial velocity

M(t): Mass as a function of time

dM/dt: Derivative of Mass with respect to time (In general, this could be a function of time.)

Ve: Nozzle velocity with respect to ground (In general, this could be a function of time.)


I gave a symbolic solution for v(t) in a previous post.

I don't know how to explain it any differently: Saying that you want to "write a differential equation in C" means absolutely nothing to me. Nothing. I also have no idea where a function called "eq.c" fits into the grand scheme of things. None of my NR version 2 or version 3 distributions have anything called eq.c (as far as I can tell).

So, I will try one more time:

If you want to do something numerically, you have to have some numbers (yes, really: numbers) and some specific equations and formulas. The equations and formulas have to have specific values for the rocket parameters.

For our first order ordinary differential equation given in the form dv/dt = f(t), and given v(0), I will try to integrate numerically from zero to values of t with the NR Simpson's rule function qsimp(). The only region of interest that I can see for the basic rocket equation is the thrust phase, since after that, the velocity will be constant.

So...

You need to feed qsimp() rule a function that can calculate dv/dt for whatever values of t you are interested in.

First of all, I would recommend making the program as generic as I could, since you may want to add or change details as you go along.

So, I would probably define a struct that holds the various parameters that define the rocket. You can test the program with hard-coded values in the struct, but eventually you might want to read the values from a file.

Here's the struct that I might define. These contain all parameters that must be defined in order to be able to get some actual numbers out of the program.


typedef struct
{
double V0; /* Initial velocity at beginning of thrust phase */
double Mf; /* Initial mass of everything = Mass of (rocket + water) */
double Mr; /* Final mass of everything = Mass of empty rocket */
double Ve; /* Nozzle velocity (exhaust with respect to the ground) */
double dMdt; /* Rate of change in mass during the thrust phase */
} Rocket;



If you are going to use qsimp, you need a function that calculates dv/dt as a function of t. For the NR2 C distribution, the function must have a float return type and have a single parameter, which is a float. In NR3, the data types are doubles, not floats.

I will mention in passing that, I just about always do all of my calculations as doubles instead of floats, and, in fact, if I were stuck with using the old versions of the Numerical Recipes routines, I would probably edit functions such as qsimp so that they would use doubles instead of floats.

Anyhow...

In order for the function to be able to use the rocket's parameters, you can have a global Rocket variable (the struct) and define the dv/dt something like:

/*
* A global variable that holds the characteristics
* of the rocket.
*/
Rocket rocket;

/* This is the basic rocket equation that we want to integrate */
float dvdt(float t)
{
return -rocket.Ve * delta_mass(t) / mass(t);
}


Now, we need to define the functions delta_mass(t) and mass(t)

Here's a possibility for mass(t):

/* Calculate mass as a function of time for constant dM/dt */
/*
* During thrust phase it is equal to
* Initial Mass - (Rate of change in Mass) * Time
* After the thrust phase, it's equal to the Mass of the rocket alone.
*/
double mass(double t)
{
double retval = rocket.Mf - rocket.dMdt * t;
if (t < 0) {
retval = rocket.Mf;
}
else if (retval < rocket.Mr) {
retval = rocket.Mr;
}
return retval;
}


And here's a possibility for delta_mass(t):

/*
Rate of change of mass with respect to time
It is constant during the thrust phase and zero afterwards.
How do we know that it is after the thrust phase?
The total mass will be equal to the mass of the rocket alone.
*/
double delta_mass(double t)
{
double retval;
if (mass(t) > rocket.Mr) { /* It's still thrusting */
retval = rocket.dMdt;
}
else { /* It's spent */
retval = 0;
}
return retval;
}


Now, the question remains: What values should we use as rocket parameters? I tried to impress upon you that we need some numbers (yes, real numbers) to actually run the program.

So, what values do you think you would use? This is important, not only just to have some numbers to plug in, but to have a feel for what the answers should be.

So, as an experiment, here's some that we can try. Maybe they are similar to what you will actually encounter with a small rocket, and maybe not, but let's try these:


rocket.V0 = 0.0; /* Initial velocity */
rocket.Mf = 1.6; /* Mass of full rocket (rocket + water) */
rocket.Mr = 0.1; /* Mass of empty rocket (rocket only) */
rocket.Ve = -1.5; /* Nozzle velocity with respect to ground */
rocket.dMdt = 6.0; /* "Burn rate" of water */



I have attached plots for dv/dt and v(t) obtained from the analytic solution.

Here's the complete program:
/*
NR2 C program showing how to use qsimp to integrate the
rocket equation
*/

#include "nr.h"
#include <math.h>

typedef struct
{
/* See the definition above. */
} Rocket;

/*
* A global variable that holds the characteristics
* of the rocket.
*/
Rocket rocket;

/* Calculate mass as a function of time for constant dM/dt */
/*
* During thrust phase it is equal to
* Initial Mass - (Rate of change in Mass) * Time
* After the thrust phase, it's equal to the Mass of the rocket alone.
*/
double mass(double t)
{
/* See the definition above. */
}

/*
Rate of change of mass with respect to time
It is constant during the thrust phase and zero afterwards.
How do we know that it is after the thrust phase?
The total mass will be equal to the mass of the rocket alone.
*/
double delta_mass(double t)
{
/* See the definition above. */
}

/* This is the basic rocket equation that we want to integrate */
/*
Note that the NR2 C program for qsimp needs floats, not doubles for
return type and argument type
*/
float dvdt(float t)
{
/* See the definition above. */
}

/*
This is the expression that comes from integrating
the rocket equation analytically.
*/
double v(double t)
{
return rocket.V0 - rocket.Ve * log(mass(0)/mass(t));
}


int main()
{
double t;
double t0;
double tf;
double deltat;
int i;
int numpoints;
/*
After development/debugging, you probably want to
change the program so that it can get the rocket
characteristics from the user command line or
from a file or whatever...
*/
rocket.V0 = 0.0; /* Initial velocity */
rocket.Mf = 1.6; /* Mass of full rocket (rocket + water) */
rocket.Mr = 0.1; /* Mass of empty rocket (rocket only) */
rocket.Ve = -1.5; /* Nozzle velocity with respect to ground */
rocket.dMdt = 6.0; /* "Burn rate" of water */

t0 = 0;
tf = 0.25;
numpoints = 100;
deltat = (tf - t0) / (numpoints - 1.0);

t = t0;
printf(" time dv(t)/dt v(t) from qsimp()\n");
for (i = 0; i < numpoints; i++) {
printf("%.5e %.5e %.5e", t, dvdt(t), v(t));
if (mass(t) > rocket.Mr)
{
double q = qsimp(dvdt,0.0, t);
printf(" %.5e", q);
}
printf("\n");
t += deltat;
}
return 0;
}


Here's the output for 10 points (not the 100 shown in the program):

time dv(t)/dt v(t) from qsimp()
0.00000e+00 5.62500e+00 0.00000e+00 0.00000e+00
2.77778e-02 6.27907e+00 1.65001e-01 1.65001e-01
5.55556e-02 7.10526e+00 3.50422e-01 3.50422e-01
8.33333e-02 8.18182e+00 5.62040e-01 5.62040e-01
1.11111e-01 9.64286e+00 8.08495e-01 8.08495e-01
1.38889e-01 1.17391e+01 1.10356e+00 1.10356e+00
1.66667e-01 1.50000e+01 1.47124e+00 1.47124e+00
1.94444e-01 2.07692e+01 1.95938e+00 1.95938e+00
2.22222e-01 3.37500e+01 2.68764e+00 2.68764e+00
2.50000e-01 9.00000e+01 4.15888e+00


A final note: There is a discontinuity on dv/dt at 0.25 seconds, since the acceleration abruptly changes to zero as the water stops being expelled. This makes the adaptive Simpson's rule not work well in case integration limits span this point. This is almost a "worst case" type of function. If we had a more complete equation (with drag and gravity), the corner might not be as sharp, and numerical integration might actually be "easier." You can write a program that changes the way it integrates as it approaches the difficult points. (I won't do it, but you can if you want to.)


Regards,

Dave

Hi Dave

I have some final questions.

How can i change the method of integration? Can i do that? E.g. the Runge Kutta 4-method and Burlisch-Stoer method.

How can i add external forces to the simulation (gravity and drag will do)?

How can i find the altitude as a function of time?

And by the way, the source for the program posted previously doesn't seem to work on my Windows-compiler. I get the following errors:

21: conflicting types for 'mass'
17: error: previous implicit declaration of 'mass' was here
33: error: conflicting types for 'delta_mass'
17: error: previous implicit declaration of 'delta_mass' was here

Thanks for your help
Best regards