dy(t)/dt = -y
In this equation t, representing time, is called the independent variable. The single unknown or dependent variable is y. The solution is strictly speaking y(t), i.e. y as a continuous function of t. Where the equation is solved numerically this will be a set of discrete values of y at corresponding values of t.
To enable o.d.e.s to be solved it is necessary to have a known fixed point in the solution space. For the present type of problem this will invariably be a starting or initial value of y at an initial time normally designated as t=0, i.e.:
y = y0 at t=0
This equation is one of the few o.d.e.s that has an analytical solution. If we solve it with initial condition:
y = 1 at t=0
y(t) = exp(-t)
To obtain a numerical solution we write:
dy/dt ~ dy/dt
Here dy and dt represent small changes in y and t rather than differentials. Hence we can calculate the change in y for a given change in t or timestep dt approximately as in general:
dy = dy/dt dt
Hence we can calculate a `new' y value since:
y(t + dt ) = y (t) + dy = y (t) + dy(t)/dt dt
For this particular equation therefore:
dy = - y dt
y(t + dt ) = y(t) - dt y(t)
Starting with a known value of y = y0 = 1
at t = 0:
y1 = y(dt) = - y0 dt + y0
y2 = y(2dt) = - y1 dt + y1
y3 = y(3dt) = - y2 dt + y2
This numerical solution is implemented here. Try changing the size of the timestep in cell B7 between say 0.01 and 0.5 and see how the numerical solution compares to the analytical one.
dy(t)/dt = - k y
dy(t)/dt = - y / T
T is called the time constant of the equation. The analytical solution of this equation is now:
y(t) = exp(-t / T )
Here is a solver for the modified problem. You can change the time constant. If you experiment with different time constants and time steps, you will observe that the accuracy of the solution depends on the ratio of the time step. Be careful at this stage to keep the time step always less than the timeconstant.
In discussion of lags in `black box' models it was shown that
the time constant of a linear o.d.e. is a measure of how
fast its solution progresses towards a final steady state.
E.g. 63% in a time equal to the time constant and 90% in
approximate 2.3 times the time constant. If we wanted to observe
99% of the response towards steady state this would require a time given
0.01 = exp (-t/T)
This gives gives a time approximately 4.6 times T.
Now make the time step just larger than the time constant, e.g. 1.1. Observe that the solution now starts to behave in a quite different, and completely incorrect oscillatory manner. If the time step is made exactly twice the time constant the solution oscillates between a maximum and minimum value. If the timestep is made any larger, then the oscillations will grow and eventually become infinite.
dx/dt = - x
dy/dt = x - y
Starting from initial values y0
and x0 we just apply
repeatedly the two relationships:
x1 = - x0 dt + x0
y1 = ( x0 - y0) dt + y0
It is clear from our investigation of stability that we must not use a time step greater than the time constant, and for reasonable accuracy it should normally be less than one tenth of the timestep. Since we would usually wish to follow most of the behaviour of the system towards its final steady state, since 4.6 times the time constant takes us 90% of the way there for a single equation a total simulation time of perhaps 5 times the time constant should be adequate,
For a single equation, or a pair of equations like those above which have the same time constant, a total of around 50 time steps would be adequate to provide sufficient accuracy and to cover the range of `interesting' behaviour.
Unfortunately, for a set of equations, the shortest time constant governs the acceptable choice of time step and longest determines the length of time for which the simulation must be run. Experiment with this using the program supplied. With the second time constant of say 5 the system will still be a long way from its final steady state, which is both variable and equal to zero, after a simulation time of 2.5. However, making the first time constant too small, i.e. less than 0.1 with the time step of 0.1, still causes instability.
This phenomenon is call stiffness and becomes a problem whenever time constants differ by more than an order of magnitude. Since in practice these may differ by several factors of 10, more sophisticated numerical methods than the one described above often have to be used to allow larger timesteps to be taken.
Next - Section 4.2: The Systematic Construction of Dynamic Models