The Mathematics of Modeling
Part 1: Time-Stepping Methods

In many ways, science is the art of predicting the future. By observing how something changes over the course of time, we learn about the conditions that affect those changes and learn to expect, and thus predict, similar changes for similar things given similar conditions.  The more information we have, the more likely it is that we can make a prediction that is correct.

Suppose I ask you to predict the weather for tomorrow.  If you've been working hard in your windowless office, and only vaguely know that it is early March, you might predict 'sunny, with some wind'.  If you had recently been outside and noticed that today's conditions were actually cloudy and windy, you would be inclined to predict much the same for tomorrow. If you were able to access an internet weather station that showed a major stormfront moving in your direction, you might intead predict 'very windy, with rain and thunderstorms'.  The conditions you know at the start of your prediction are called initial conditions.

The accuracy of a prediction will also depend on the length of time that you are trying to predict over. Continuing with our weather example - you would probably be able to make a very good prediction of the weather five minutes from now, even with only limited personal observations. Predicting twenty-four hours from now is much less likely to be accurate, and predicting a month from now becomes little more than a guess, even if you have the complete array of technology to determine your initial conditions.

Time Stepping
	Time-stepping methods try to take advantage of the idea that changes in a system can be predicted pretty well, if the prediction time is small. Presume for example that we are tracking a hurricane in the Atlantic, and that we expect landfall in twenty-four hours, but would like to predict where the hurricane will be at that point. We start by dividing our prediction time into many smaller pieces - small enough so that we think we can make a good prediction, but not too small! - the more pieces we have the more work we have to do, and it doesn't help if it takes more than 5 minutes to figure out what happens (happened!) in those five minutes! A good sized piece for this example might be an hour.

This small piece of time has a special name - it is known as a time step, and it represents a change in time.  We represent change mathematically using the greek letter delta, capitalized, so if we use t to represent time, our time step would be represented as Δt. We also label our start time with a zero or 'naught' so that 't-naught' (t\(TraditionalForm\`\_0\)) is the start time, and after the first time step we are at t\(\_1\),then t\(TraditionalForm\`\_2\), t\(TraditionalForm\`\(\(\_3\)\(,\)\)\)t\(TraditionalForm\`\(\(\_4\)\(,\)\)\)etc...  To talk about 'some time step' in general, we label the 'time now' as t\(TraditionalForm\`\_n\)("t sub n"), which will represent the beginning of a time step, and the time at 'now plus one time step', or the end of the time step, becomes t\(TraditionalForm\`\(\ \_\(n + 1\)\)\). ("t sub n plus one").

In general, time will be what we call and independent variable, meaning that we don't expect anything to affect how time goes by.  Other information we are interested in, like the position of the hurricane, might require the current time, as well as other information, to determine a value. Variables whose values depend on other variables are called dependent variables. If our variables are dependent on time, we will often label them in the same way that we label our times. If we keep track of position using x and y as variables, the initial position is the position that corresponds with the start time, and so we label those variables with a naught as well, giving starting position {x\(TraditionalForm\`\_0\), y\(TraditionalForm\`\_0\)}. The position at the beginning of any time step would then be {x\(TraditionalForm\`\_n\), y\(TraditionalForm\`\_n\)}, and at the end would be {x\(TraditionalForm\`\_\(n + 1\)\), y\(TraditionalForm\`\_\(n + 1\)\)}.

The key piece in any time stepping method is determining how things change in time. Velocity is the measurement of how position changes in time, in calculus, 
				v  =  (dx/dt),
Where v is the velocity, x is a position variable, t is time, and the d is another represention of 'the change in' (like the delta!), although with the calculus assumption that this change is infinitesimally small and that all we care about is the relationship between the two changes (in x and in t).  We are not dealing with infinitesimally small pieces, instead we are insterested in an 'average velocity', with bigger (Δ) changes:
				v  =  Δ(x/Δ)t.
So velocity becomes the 'change in position' over the 'change in time' (which is our time step).  We can rewrite the change in position as 'the position at the end of a time step' (x\(TraditionalForm\`\_\(n + 1\)\)) minus 'the position at the beginning of the time step (x\(TraditionalForm\`\_n\)). Now the velocity description looks like:
				\(TraditionalForm\`v\  = \ \(x\_n - x\_\(n + 1\)\)\/Δt\).
We assume that we know the position at the beginning and that we want to calculate the position at the end, so we manipulate the equation above to get
				(Subscript[x, n + 1]) = 
        (Subscript[x, n]) + Δ t  v .
We also know the time step (we choose it!), so the only question is what the velocity is.  Officially it is an 'average' velocity over the time, and the velocity could change quite a bit during that time!  However, we will use a method called Euler's method and assume that we know (or will have calculated) the velocity at the beginning of the time step (v\(TraditionalForm\`\_n\)), and that this velocity will be a good estimate of the average velocity for a short period of time. So now we have a time-stepping equation:
				(Subscript[x, n + 1]) = 
        (Subscript[x, n]) + Δ t  (Subscript[v, n]) .
	If we have a second position variable (like y), we would need a second velocity variable (let's use w) and a second equation:
				(Subscript[y, n + 1]) = 
        (Subscript[y, n]) + Δ t  (Subscript[w, n]) .
	How do we find the velocities? The change in velocity is knows as acceleratious, 
				a  =  (dv/dt),
so, following the same process, we can get time stepping equations for velocities, too!
			         (Subscript[v, n + 1]) = 
        (Subscript[v, n]) + Δ t  (Subscript[a, n]) ,
			(Subscript[w, n + 1]) = 
        (Subscript[w, n]) + Δ t  b\(TraditionalForm\`\(\(\_n\)\(\ \)\)\).
	Eventually we need to stop making new variables, and in position problems, that usually occurs with the acceleration.  If we were dealing with gravity, we might say that the accelerations are constant - for every time step there might be no acceleration in the x direction and acceleration equal to gravity in the negative (downwards) y direction.
				(Subscript[a, n])  =  0
(Subscript[b, n])  =  -9.81
	One thing to notice is that for every variable we have, we will need an initial condition! We already know that we need starting positions, (Subscript[x, 0])and (Subscript[y, 0]), but we will also need starting velocities  (Subscript[v, 0])and (Subscript[w, 0]), as well as starting accelerations (although we may be able to calculate those if we have formulas). 
	After we have calculated the values at the end of each time step, we use those values for the beginning of the next time step, so what was (Subscript[x, n + 1])is stored as (Subscript[x, n]),(Subscript[v, n + 1])becomes (Subscript[v, n])and so forth.Let's look at an example...

Projectile Motion
	We can do a simple time-stepping simulation of the motion of a shotput. First we define our initial positions, in meters, presuming that the shotputter is average height.

(Subscript[x, 0]) =  0;
(Subscript[y, 0]) =  2;

Then we estimate initial values for velocity (in meters per second) - we will make them equal assuming a 45 degree release.

(Subscript[v, 0]) =  4;
(Subscript[w, 0]) =  4;

Our time step will be one tenth of a second, and our accelerations will be constant.

Δt  =  0.1;
a  =  0;  b  =  -9.81;

Now we make Euler method function:

EulerCalc[val_, rate_, timestep_]  :=  
        val  +  timestep*rate;

We then set our n values to the starting values:

(Subscript[x, n]) =  (Subscript[x, 0]); (Subscript[y, n]) =  (Subscript[y, 0]); (Subscript[v, n]) =  (Subscript[v, 0]); (Subscript[w, n]) =  (Subscript[w, 0]); 
    (Subscript[t, n]) =  0;

To calculate the first step, we use Euler's method for each variable, and add a time step to the time.

(Subscript[x, n + 1])  =  
        EulerCalc[(Subscript[x, n]), (Subscript[v, n]), Δt];
(Subscript[y, n + 1])  =  
        EulerCalc[(Subscript[y, n]), (Subscript[w, n]), Δt];
(Subscript[v, n + 1])  =  
        EulerCalc[(Subscript[v, n]), a, Δt];
(Subscript[w, n + 1])  =  
        EulerCalc[(Subscript[w, n]), b, Δt];
(Subscript[t, n + 1]) =  (Subscript[t, n]) +  Δt;

Then we output the position results

{(Subscript[x, n + 1]), (Subscript[y, n + 1])}

We then take all the values at n+1 (the 'end' of the first time step) and store them in the n variables (the 'beginning' of the next time step).

(Subscript[x, n]) =  (Subscript[x, n + 1]); (Subscript[y, n]) =  (Subscript[y, n + 1]); (Subscript[v, n]) =  (Subscript[v, n + 1]); 
    (Subscript[w, n]) =  (Subscript[w, n + 1]); (Subscript[t, n]) =  (Subscript[t, n + 1]);

And repeat the process

(Subscript[x, n + 1])  =  
        EulerCalc[(Subscript[x, n]), (Subscript[v, n]), Δt];
(Subscript[y, n + 1])  =  
        EulerCalc[(Subscript[y, n]), (Subscript[w, n]), Δt];
(Subscript[v, n + 1])  =  
        EulerCalc[(Subscript[v, n]), a, Δt];
(Subscript[w, n + 1])  =  
        EulerCalc[(Subscript[w, n]), b, Δt];
(Subscript[t, n + 1]) =  (Subscript[t, n]) +  Δt;

{(Subscript[x, n + 1]), (Subscript[y, n + 1])}

This is obviously a good use for repeating and storing the result using the Table function.

(Subscript[x, n]) =  (Subscript[x, 0]); (Subscript[y, n]) =  (Subscript[y, 0]); (Subscript[v, n]) =  (Subscript[v, 0]); (Subscript[w, n]) =  (Subscript[w, 0]); 
    (Subscript[t, n]) =  0;
Trajectory  =  
      Table[
((Subscript[x, n + 1])  =  
            EulerCalc[(Subscript[x, n]), 
              (Subscript[v, n]), Δt]; 
(Subscript[y, n + 1])  =  
            EulerCalc[(Subscript[y, n]), 
              (Subscript[w, n]), Δt]; 
(Subscript[v, n + 1])  =  
            EulerCalc[(Subscript[v, n]), 
              a, Δt]; 
(Subscript[w, n + 1])  =  
            EulerCalc[(Subscript[w, n]), 
              b, Δt]; 
(Subscript[t, n + 1]) =  
            (Subscript[t, n]) +  Δt; 
(Subscript[x, n]) =  
            (Subscript[x, n + 1]); (Subscript[y, n]) =  (Subscript[y, n + 1]); (Subscript[v, n]) =  (Subscript[v, n + 1]); 
          (Subscript[w, n]) =  (Subscript[w, n + 1]); 
          (Subscript[t, n]) =  (Subscript[t, n + 1]); 
{(Subscript[x, n + 1]), 
            (Subscript[y, n + 1])}), {i, 12}]

Which can then be plotted...

ListPlot[Trajectory, Prolog -> AbsolutePointSize[3]]

The Stella software is particularly good for time stepping routines, as it determines the formulas from the stocks and flows, using the converters to make the relationships.  You can see these relationships if you look at the equation (bottom) level of a Stella model, where you will see equations that include a time step, labelled dt in Stella. The key behind accurate modeling is ensuring that you have these relationships correct - if you do, you too, can 'predict the future'!

This notebook created by Dr. Garrett R. Love for the Shodor Education Foundation. Copyright 2001