Notebook 1

Introduction to Simulation and Modelling

Acknowledgements

Parts of the material in this section was taken from a number of sources, which are listed below

Mathematica 4.0, S. Wolfram.
Mathematical Modelling and Mathematica, J.E. Cruthirds and L.R. Hitt, Department of Mathematics and Statistics, University of South Alabama, material taken from the Web.
An introduction to Mathematica, N. Blachman, material taken from the Web.

The introductory text on modelling and simulation was taken from: 
Sloot, P.M.A., "Modelling and Simulation," in Proceedings of the 1994 CERN School of Computing, CERN, ISBN 92-9038-069-7 CERN press, ISSN 0304-2898: Also Technical Report: CS-94-12, 1994.

See also reference at the end of Chapter I.

Introduction

This notebook will  introduce some concepts of modelling and simulation.   **It  is not  required to hand this  notebook  in.**

Important  instructions  for  all students:
Each notebook  contains several questions that must be handed in. Theses questions are labeled accordingly:

Required: These questions must be done by all students.

Advanced:  These  are considered to be more difficult and may be assigned  by the instructor.

Optional:   These questions are optional  to all students  and will be counted towards   extra  marks.

In addition to studying this material, you are supposed to look at the notebook Appendix.nb where an introduction to Mathematica is given. At the end of that notebook  you will  find a number of exercises.

I Introduction to Simulation and Modelling

I.1 The Computer  Experiment

The use of computers as an "experimentation environment" is becoming increasingly important in a wide area of applications. On the one hand, computers are used in experimental set-ups  such as a measuring, controlling, and data analysis devices, which are essential  to the accurate measurement and data handling. On the other  hand, we have the field where computers are used to perform some simulation of a "real-world"  phenomenon. One of the crucial components of that research field is the correct abstraction of a real-world  phenomenon to a conceptual model and the translation  into a computational  model that can be validated. This leads us to the notion of a computer experiment, where the model and the computer take the place of the "classical" experimental set-up where the simulation replaces the experiment.
In general, computer experiments are done for the following   purposes:

I.2 Systems, Models, and Simulation

I.2.1 Systems and System Environment

To model a system, it is necessary to understand the concept of a system and the system constraints. We define a system as a group of objects that are joined together in some interaction or interdependence toward the accomplishment of some purpose. An example is a production system manufacturing automobiles. The machines, component parts, and workers operate jointly along an assembly line to produce a high-quality    vehicle.
A system is often affected by changes occurring outside the system. Such changes are said to occur in the system environment. In modelling systems, it is necessary to decide on the boundary between the system and its environment. This decision may depend on the purpose of the study.

I.2.2 Components of a System

In order to understand and analyse a system, we need to define a number of terms. 
i) entity:  This is an object of interest  in the system. An 	     attribute  is a property of an entity. 
ii) activity:  Represents  a time period of specified length. If a bank is being studied, customers might be one of the entities, the balance in their checking accounts might be an attribute, and making deposits might be an activity.
iii) state of a system:  A collection of variables necessary to describe a system at a particular time, relative to the objectives of a study. In a study of a bank, examples of possible state variables are the number of busy tellers, the number of customers in the bank, and the time of arrival of each customer in the bank.
We categorise systems to be one of two types, discrete and continuous. 

iv) Discrete system: The state variables change instantaneously at separated points in time. A bank is an example of a discrete system, since state variables, e.g. the number of customers in the bank, change only when a customer arrives or when a customer finishes being served and departs. 
v)  Continuous system: The state variables change continuously with respect to time. An aeroplane moving through the air is an example of a continuous system, since state variables such as position and velocity can change continuously with respect to time. 

Few systems in practice are completely discrete or completely continuous, but since one type of change predominates for most systems, it will usually be possible to classify a system as being either discrete or continuous.
With  respect to discrete systems, we define an event as an instantaneous occurrence that may change the state of the system. The term endogenous is used to describe activities and events occurring within a system, and the term exogenous is used to describe activities and events in the environment that  affect the system.

I.2.3 Model  of a System

A model is defined as a representation of a system for the purpose of studying the system. In practice, what is meant by "the system" depends on the objectives of a particular study. For most studies, it is not necessary to consider all the details of a system; thus, a model is not only a substitute for a system, it is also a simplification of the system. However, there should be sufficient detail in the  model to permit valid conclusions to be drawn about the real system.
Different  models of the same system may be required as the purpose of investigation  can change. For example,  one may want to study a bank to determine the number of tellers needed to provide adequate service for customers who want just to cash a check or make a savings deposit. The model can be defined to be that portion of the bank consisting of the tellers and the customers waiting in line or being served. If the loan officer and the safety deposit boxes are to be included, the definition of the model must be expanded in an obvious way.
Just as the components of a system were entities, attributes, and activities, models are represented similarly. However, the model contains only those components that are considered to be relevant to the study.

Table 1.1: Overview of simulation concepts.

I.2.4 Experimentation and Simulation

At  some point in the lives of most  systems, there is a need to study them to try to gain some insight into the relationship  among various components, or to predict performance under some new conditions being considered. Figure 1.1 maps out different  ways in which a system might be studied.

Figure 1.1: Ways to study a system.

I.2.4.1 Experiment with the Actual System  vs. Experiment with a Model of the System

If it is possible (and cost-effective)  to alter the system physically and then let it operate under the new conditions, it is probably desirable to do so. In this case there is no question about whether what we study is relevant.  However, it is rarely  feasible to do this, because such an experiment would often be too costly or too disruptive to the system. For example, a bank may be contemplating  reducing the number of tellers  to decrease costs, but actually trying this could lead to long customer delays and queues. In more abstract terms this "system" might not even exist, but we nevertheless want to study it in its various proposed alternative  configurations to see how it should be built in the first place. Examples of this situation might be modern flexible manufacturing facilities, or strategic nuclear weapons systems. For these reasons, it is usually necessary to build a model as a representation  of the system and study it as a surrogate for the actual system. When using a model, there is always the question of whether it accurately reflects the system for the purposes of the decisions to be made.

I.2.4.2 Physical  Model vs. Mathematical  Model

To most people, the word "model" evokes images of clay cars in wind tunnels, cockpits disconnected from their aeroplanes to be used in pilot training, or miniature super tankers scurrying about in a swimming pool. These are examples of physical models (also called iconic models), and are not typical of the kinds of models that  are usually of interest in systems analysis and computer simulations. In some circumstances it has been found useful to build physical models to study engineering or management systems; examples include tabletop scale models of material-handling systems, and in at least one case a full-scale physical model of a fast-food restaurant inside a warehouse, complete with full-scale, real (and presumably hungry) humans. 

The vast majority of models built for such purposes are mathematical, representing a system in terms of logical and quantitative relationships that are then manipulated and changed to see how the model reacts, and thus how the system would react—if the mathematical model is a valid one. 

Perhaps the simplest example of a mathematical model is the familiar  relation  d = vt, where v is the speed of travel, t is the time spent travelling, and d is the distance travelled. This might provide a valid model in one instance (e.g., a space probe to another planet after it has attained its flight velocity) but a very poor model for other purposes (e.g., rush-hour commuting on congested urban freeways).

I.2.4.3 Analytical Solution vs. Simulation

Once we have built a mathematical model, it must then be examined to see how it can be used to answer the questions of interest about the system it is supposed to represent. If the model is simple enough, it may be possible to work with its relationships and quantities to get an exact, analytical solution. In the d = vt example, if we know the distance to be travelled and the velocity, then we can work with the model to get t = d/v as the time that will be required. This is a very simple, closed-form solution obtainable with just paper and pencil, but some analytical solutions can become extraordinarily complex, requiring vast computing resources; inverting a large non sparse matrix is a well-known example of a situation in which there is an analytical formula known in principle, but obtaining it numerically in a given instance is far from trivial. If an analytical solution to a mathematical model is available and is computationally efficient, it is usually desirable to study the model in this way rather than via a simulation. However, many systems are highly complex, so that valid mathematical models of them are themselves complex, precluding any possibility of an analytical solution.  In this case, the model must be studied by means of simulation, i.e., numerically exercising the model for the inputs in question to see how they affect the output measures of performance.
While there may be an element of truth to pejorative old saws such as "method of last resort" sometimes used to describe simulation, the fact is that we are very quickly led to simulation in many situations, due to the sheer complexity of the systems of interest and of the models necessary to represent them in a valid way.
Given, then, that we have a mathematical model to be studied by means of simulation (henceforth referred to as a simulation model), we must then look for particular tools to execute this model (i.e. actual simulation). It is useful for this purpose to classify simulation models along three different  dimensions:

I.2.4.4  Static vs. Dynamic Simulation Models

A static simulation  model is a representation of a system at a particular  time, or one that may be used to represent a system in which time simply plays no role; examples of static simulations are Monte Carlo models (see chapters 2 and 7). On the other hand, a dynamic simulation  model represents a system as it evolves over time, such as a conveyor system in a factory.

I.2.4.5 Deterministic  vs. Stochastic Simulation Models

If a simulation  model does not contain any probabilistic (i.e., random) components, it is called deterministic. A complicated (and analytically intractable) system of differential equations describing a fluid flow is an example of such a model. In deterministic models, the output is "determined" once the set of input quantities and relationships in the model have been specified, even though it might take a lot of computer time to evaluate what it is. 
Many systems, however, must be modelled as having at least some random input components, and these give rise to stochastic simulation models. Most queuing and inventory systems are modelled stochastic. Stochastic simulation models produce output that  is itself  random, and must therefore be treated as only an estimate of the true characteristics of the model; this  is one of the main disadvantages of simulation.

I.2.4.6 Continuous vs. Discrete Simulation Models

Loosely speaking, we define discrete and continuous simulation  models analogously to the way discrete and continuous systems were defined above. It should be mentioned that a discrete model is not always used to model a discrete system and vice versa. The decision whether to use a discrete or a continuous model for a particular system depends on the specific objectives of the study. For example, a model of traffic flow on a freeway would be discrete if the characteristics and movement of individual cars are important.    Alternatively, if the cars can be treated "in the aggregate," the flow of traffic can be described by differential equations in a continuous  model.

I.2.5 A Closer Look at System Models

Although many attempts have been made throughout the years to categorise systems and models, no consensus has been arrived at. However, it is convenient to make the following distinction between the different models:

I.2.5.1 Continuous-Time Models

Here the state of a system changes continuously over time. These types of models are usually represented by sets of differential equations. A further subdivision would be:
* Lumped parameter models expressed in ordinary differential equations (ODE's):

* Distributed parameter models expressed in partial differential equations (PDE's):

Figure 1.2: Trajectory of continuous-time model.

I.2.5.2 Discrete-Time Models

With discrete-time models, the time axis is discretised. The system state changes are commonly represented by difference equations. These types of models are typical to engineering systems and computer-controlled systems. They can also arise from discrete versions of continuous-time models.

The time-step used in the discrete-time model is constant.

Figure 1.3: Trajectory of a discrete-time model.

I.2.5.3 Discrete-Event Models

In discrete-event models, the state is discretised and "jumps" in time. Events can happen any time but only every now and then at (stochastic) time intervals. Typical examples come from "event tracing" experiments, queuing models, Ising spin simulations, image restoration, combat simulation, etc.

Figure 1.4: Trajectory of discrete-event model.

I.3 Modelling and Simulation Cycle

During the development of a computer experiment, we can identify three major phases, each with its own challenges:
Modelling  Phase
The first step to simulation is the development of an abstract model of the real-world  system under study. Strange enough there are hardly any formal methods of modelling which supports completeness, correctness, and efficiency. The lack of such methods might harm the success of this research field significantly. There are complete journals and internet sites dedicated to the compilation of the consequences of using insufficient, incorrect and inefficient models. For instance some flight simulation experiments were reported where the underlying simulator worked perfectly, but where deficiencies in the CAD/CAM model resulted in catastrophic misinterpretation of the simulations. Some promising new approaches in this field are the use of structural interactive modelling methods (such as presented in STELLA) and the concept of bond graph modelling used in control systems design or object-oriented modelling.
Simulation Phase
Solvers form the kernel of simulation systems. Here we refer to mathematical methods that make the underlying real-world  models discrete. A rough distinction can be made between solvers for discrete (event) systems and solvers for continuous systems. Solvers for discrete event systems are event-set algorithms, which ensure that the events occur in the proper order and at the proper time. The model system is tracked in time. The conventional solvers for continuous systems are finite difference, finite element/volume, and a large class of linear algebra solvers. Also of special interest are the new results obtained with particle methods.
Figure 1.5 shows  how various phases in modelling and simulation  are bound together.

Figure 1.5: Overview of conceptual steps in simulation and modelling cycle.

Computational Phase
In this phase we concentrate on the mapping of the different solvers to the machine architecture. Since the type of problems we are interested in are computationally very demanding, much research effort is going on in the efficient use of modern architectures (parallel computer systems) for simulation. As an example, Fig. 1.6 shows the memory requirements (in Words) and computational performance (in Floating-Point Operations per Second) needed to investigate a number of challenging problems by simulation.

Figure 1.6: Examples of challenging computational real-world research topics and their computational requirements.

In Fig. 1.7 we summarise the various conceptual steps required to map a real-world system onto a computer.

Figure 1.7: Summary of the various conceptual stages in the development of a computer experiment.

The computation specific model consists of representing the derived conceptual model into a language that can be implemented. One of the things that often goes wrong in these stages is that the scientist tends to confuse his code with a conceptual model or even worse, with the real-world system itself. This is a situation we should be prepared for and try to avoid at all costs. One way to maintain a critical attitude is to carefully design and test sub stages in the modelling and simulation cycle and to add primitives to the simulation that constantly monitor for "unrealistic" results.

I.4 Model Execution: Event driven versus Time driven

We have seen that in continuous systems the state variables change continuously with respect to time, whereas in discrete systems the state variables change instantaneously  at separate points in time. Unfortunately for the computational experimentation there are but a few systems that are either completely discrete or completely continuous state, although often one type dominates the other in such hybrid systems. The challenge here is to find a computational model that mimics closely the behaviour of the system, specifically the simulation time-advance approach is critical.
If we take a closer look into the dynamic nature of simulation models, keeping track of the simulation time as the simulation proceeds, we can distinguish between two time-advance approaches: time-driven and event-driven.

Time-driven   Simulation
For continuous systems, time-driven  simulations  advance time with a fixed increment. With this approach the simulation  clock is advanced in increments of exactly Δt time units. Then after each update of the clock, the state variables are updated for the time interval [t, t+Δt]. This is the most widely known approach in simulation of natural systems. 
Less widely used is the time-driven  paradigm applied to discrete systems. In this case we have specifically to consider whether:
• The time step Δt is small enough to capture every event in the discrete system. This might imply that we need to make Δt arbitrarily small, which is certainly not acceptable with respect to the computational times involved.
• The precision required can be obtained more efficiently through the event-driven execution mechanism. This primarily means that we have to trade efficiency for precision.

Event-driven  Simulation
In event-driven  simulation the next-event time advance approach is used. For  the case of discrete systems this  method consists of the following  phases:

The advantage of this approach is that periods of inactivity  can be skipped over by jumping the clock from event time to the next event time. This is perfectly safe since by definition all state changes only occur at event times. Therefore causality is guaranteed. The event-driven approach to discrete systems is usually exploited in queuing and optimisation problems. However, as we will  see next, it is often also a very interesting paradigm for the simulation of continuous systems. 
Consider a continuous system where every now and then (possibly at irregular or probabilistic time steps) discontinuities  occur, for instance in the temperature of a room where the heating is regulated  in some feed-back loop:

Figure 1.8: Typical discontinuities  in Time versus State trajectories of continuous systems or its (higher order) derivative with respect to time.

The difficulty of an efficient time-driven simulation of such a system is in the  integration method applied. Specifically, multi-step  integration  methods to solve the underlying differential equations might prove not to work in this case. The reason is that these methods use extrapolation to estimate the next time step, this would mean that they will try to adapt to the sudden change in state of the system thus reducing the time step to infinite small sizes.
Other examples are:
Continuous Ising spin systems: A configuration is defined by the spin variables s(ν) = ± 1, specified at the vertices ν of a two or three dimensional lattice and an independent Poisson process of attempted spin change arrivals are associated with the vertices of the lattice. If in such a system an attempted  change arrives at vertex ν , the spin s(ν) is changed to -s(ν) with probability P, and with probability 1-P the spins remain unchanged.
Billiard ball models: Used in gas dynamics, where N equal balls  move without  friction in 2D and where each ball moves along a straight line with constant velocity, until  it comes into contact with  another ball. Upon such a contact a collision  occurs where upon the two participating balls instantly  change their velocities and directions. 
It is clear that time-driven  simulation might not prove to be acceptable since this mechanism is slow and/or not sufficiently precise if Δt is not chosen to be very small. An event-driven  simulation where the state of the balls is updated from collision to collision of any ball, might prove much more accurate and efficient.
Let's explore this example in a little more detail.

I.4.1 Billiard Balls - Example

Consider a Billiard Ball model to simulate the gas dynamics of 4 gas molecules (as long as we do not try to do any physics with this simple system it can perfectly well explain the concepts we are discussing). Assume that the molecules move along the straight line with constant velocity until they hit the wall of a 2D system (kinetic energy is conserved in this model). We can now investigate the consequences of the two paradigms time-driven versus event-driven simulation:

input stop simulation time, Δt;
initialise balls with x,y position and velocity vx,vy in x,y direction;
time = 0.0;
while (time < stop simulation)
{
for each ball do {
		x += Δt * vx;
		y += Δt * vy;
		if ((x == 0)||(x == 1)) vx = -vx;
		if ((y == 0)||(y == 1)) vy = -vy;
	}
time += Δt;
}

Code 1.2: Pseudo code for time-driven gas simulation.

For  a simulation  time of 1.0 (arbitrary  units) with  a Δt  of 0.001, implying 1000 state evaluations, this time-driven approach (see Table 1.2) results in the following simulated  behaviour of the 4 balls  over time:

Figure 1.9: Time-driven state evaluation of gas model, 1000 state evaluations, 4 molecules.

Next we consider the case of an event-driven simulation of this 4-molecule  gas. The pseudo code shown in Table 1.3 takes into account all the typical aspects of detecting, handling and processing events within the correct time frame.

input stop_simulation_time;
initialise balls with x,y position and velocity vx,vy in x,y direction;
for each ball do {
	impact_time = collision_time(ball);
	schedule(MOVE, impact_time, ball);
}
prev_time = 0.0;
while (time() < stop_simulation_time) {
	next_event(&event, &ball);
	switch(event) {
	case MOVE:
		update_positions(time() - prev_time);
		impact_time = collision_time(ball);
		schedule(MOVE, impact_time, ball);
		prev_time = time();
		break;
	}
}
collision_time(ball)
{
	if (vx > 0) {
		t0 = (1 - x) / vx;
	} 
	else {
		t0 = -x / xv;	
	}
	if (vy > 0) {
		t1 = (1 - y) / vy;
	}
	else {
		t1 = -y / vy;	
	}
	return min(t0, t1);
}
update_positions(Δt)
{
	for each ball do {
		x += Δt * vx;
		y += Δt * vy;
		if (x == 0 || x == 1)
			vx = -vx;
		if (y == 0 || x == 1)
			vy = -vy;
	}
}

Code 1.3: Pseudo code for an event-driven gas simulation.

For a simulation time of 1.0 (arbitrarily units) with 60 state evaluations, this event-driven approach results in the following simulated behaviour of the 4 balls over time:

Figure 1.10: The trajectory of an event driven billiard ball simulation
From these examples we see that the cumulated error over time in the time-driven simulation, even for a Δt as small as 0.001, results in a violation of the conservation of momentum. The gain we have in event-driven simulations is however partly frustrated by the growth in the (complexity of the) simulation code. If we would like to go to modern (parallel) architecture's, event-driven simulation  becomes increasingly more difficult to code due to the complicated synchronisation protocols. There is always the necessity to find a balance between the pro's and con's of time-driven versus event-driven simulation.
In general we observe that:

•	If the time between the succeeding events becomes small, time driven is preferred (in the example this situation would occur for a system of billions of molecules to be simulated). This is the frequency parameter.
•	The overhead per event (i.e., the state update) is larger with event-driven simulation. This is referred to as the event overhead.
•	The amount of work between any two events plays an important role. This we call the granularity of the system.

The factors: frequency, overhead and granularity, together with the consequences for implementation on modern architecture's and the required accuracy, must be studied carefully before making a decision on what simulation mechanism should be used.

I.5 Analysis of Simulation Results

During the development of the simulation model, we must ensure that the model is correctly implemented and that it is representative of the real system. These two steps are called model verification and model validation, respectively. After the model development is complete, the next two issues we will face are those of deciding how many of the initial observations should be discarded to ensure that the model has reached a steady state and how long to run the simulation. These issues are called transient removal and stopping criterion, respectively. These four issues are the main topics of this section.

I.5.1 Model Verification Techniques

The accuracy of a simulation model is measured by the closeness of the model output to that of the real systems. Since a number of assumptions about the behaviour of real systems are made in developing the model, there are two steps in measuring the accuracy. The first step is checking whether or not the assumptions are reasonable. The second step is whether or not the model implements those assumptions correctly. These two steps are called validation and verification, respectively. Validation is concerned with the representativeness of the assumptions, and verification is related to the correctness of the implementation. Verification can also be called debugging, that is, ensuring that the model does what it is intended to do.
In the programming language literature, a number of techniques can be found for debugging. Any combination of these techniques can be used to verify the model.

I.5.2 Model Validation Techniques

Validation  refers to ensuring that the assumptions used in developing the model are reasonable in that, if correctly implemented, the model would produce results close to that observed in real systems. The validation techniques depend upon the assumptions and, hence, on the system being modelled. Thus, unlike verification techniques that are generally applicable, the validation techniques used in one simulation may not apply to another.
One of the classical papers in the validation literature  is that of Naylor and Finger (1967), where a three-step approach is given for validating  a simulation  model:

The primary objective during the first step is to develop a model with high face validity, i.e., a model that, on the surface, seems reasonable to people who are knowledgeable about the system under study. In order to develop such a model, the simulation modellers should make use of all existing information, including the following:

Test the Assumptions of the Model Empirically.
The goal of the second step of validation is to test quantitatively the assumptions made during the initial stages of model development.
One of the most useful tools during the second step of validation is sensitivity analysis. This can be used to determine if the simulation output changes significantly when the value of an input parameter is changed, when an input probability distribution is changed, or when the level of detail for a subsystem is changed. If the output is sensitive to some aspect of the model, then that aspect must be modelled  carefully.
Determine How Representative the Simulation Output Data Are
The most definitive test for a simulation  model's validity is establishing that its output data closely resembles the output data that would be expected from the actual (proposed) system. If a system similar to the proposed one exists, then a simulation model of the existing system is developed and its output data are compared to those from the existing system itself. If the two sets of data compare "favourably", then the model of the existing system is considered "valid". The model is then modified so that it represents the proposed system. The greater the communality between the existing and proposed systems, the greater our confidence in the model of the proposed system. There is no completely definitive approach for validating the model of the proposed system. If there were, there might be no need for a simulation  model in the first place.
If there is no an existing system similar to the proposed system or if there is an existing system but no definitive output data, then it is still worthwhile to have system experts review the simulation output data for reasonableness. (Care must be taken in performing this exercise, since if one knew exactly what output to expect, there would be no need for a model.) Animation may also be an effective way for experts to evaluate the validity of a simulation  model.
Up to now we have discussed validating a simulation model relative to past or present system output data; however, a perhaps more definitive test of a model is to establish its ability to predict future system behaviour, when the model input data match the real inputs and when a policy implemented in the model is implemented at some point in the system. In other words, the structure of the model should be accurate enough for the model to make good predictions, not just for one input data set, but for the range of input data sets, which are of interest.

I.5.3 Transient Removal

In most simulations, only the steady-state  performance, that is, the performance after the system has reached a stable state, is of interest. In such cases, results of the initial part of the simulation should not be included in the final computations. This initial part is also called the transient state. The problem of identifying the end of the transient state is called transient removal.
The main difficulty with transient removal is that it is not possible to define exactly what constitutes the transient state and when that transient state ends. All methods for transient state removal are therefore heuristic. Examples of such methods are: (i) simply use very long runs; that is, runs that are long enough to ensure that the presence of initial conditions will not affect the result; or (ii) use proper initialisation, which requires starting the simulation in a state close to the expected steady state; and (iii) various other methods that are based on the assumption that the variability during the steady state is less than during the transient state, which is generally true.

I.5.4 Stopping Criteria: Variance Estimation

It is important that the length of the simulation be properly chosen. If the simulation is too short, the results may be highly variable. On the other hand, if the simulation is too long, computing resources and manpower may be unnecessarily wasted. The simulation should run until the confidence interval for the mean response narrows to a desired width.

I.6 References

Fishwick, P.A., Building Digital Worlds (Prentice-Hall, 1995).
Gaylord, R.J. and Wellin, P.R., Computer Simulations with Mathematica (Springer-Verlag, 1994).
Gould, H. and Tobochnik, J., Computer Simulation Methods part I and II (Addison Wesley, 1987).
Koonin, S.E. and Meredith, D.C., Computational Physics (Addison Wesley, 1989).
Naylor, T.H. and Finger, J.M., "Verification of Computer Simulation Models," Management Science 14, 92-101 (1967).
Overeinder, B., Hertzberger, L.O., and Sloot, P.M.A., "Parallel Discrete Event Simulation," in Proceedings of the Third Workshop Computersystems, Ed. Withagen, W.J., Faculty of El. Eng., Eindhoven University, The Netherlands, 1991, pp. 19-30.
Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T., Numerical recipes in C (Cambridge University Press, Cambridge, 1988).
Ross, S.M., A course in Simulation  (Maxwell MacmillianNew York, 1991).
Sloot, P.M.A., "Modelling and Simulation," in Proceedings of the 1994 CERN School of Computing, CERN, ISBN 92-9038-069-7 CERN press, ISSN 0304-2898: Also Technical Report: CS-94-12, 1994.
Sloot, P.M.A., Parallel Scientific Computing and Simulation, Tech. Rept. CS-95-xx , January, 1995.
Sloot, P.M.A., Computational Physics, Technical Report CS95-06
Weiss, G.H., "Random numbers and their Applications," American Scientist 71, 65-71 (1983).
Wolfram, S., Mathematica: A system for doing mathematics by computer (Addison Wesley, 1991).