Module 4: Simulating Continuous Systems

Simulating a simple chemical reaction

We'll start by modeling a simple reaction involving three molecules:

One molecule of substance A reacts with one molecule of substance B to produce one molecule of C.
Similarly, a molecule of C can breakdown and produce one molecule each of A and B.
We will write this as two reactions:
1. A + B → C, and
2. C → A + B

Exercise 1: Download and unpack this ZIP file into a new directory. Upon unpacking, you will see that it creates a new subdirectory called VirtualReactor. Go into that subdirectory, then compile and execute the VirtualReactor.java file. This will bring up a window. Here's what you do:

Note that the Model menu offers you the choice of different models. Select the first one: standard simulation.
In the section entitled "Initial numbers of molecules" enter 100 in the box for A, 50 in the box for B, and 0 for C. Then click "Change" in the same row.
Hit the Reset button in the "Animation" section.
Then hit the "Go" button, and simulate until the time shows about 100 (this is in the rightmost label of the "Stats" section).
About time 100, hit "Stop" and then "Plot". You should see the time-evolution of the concentrations of each of A, B and C.
Finally, compare how long it takes to simulate 100 time units when there are large numbers of molecules (10000 of A, 5000 of B, 0 of C).

Now let's explain how the software works:

It simulates molecules in two dimensions, with small blue "A" molecules, slightly larger green "B" molecules, and much larger yellow "C" molecules.
The molecules move about randomly, simulating molecules in a solution or in a gaseous state.
At any given instant, exactly one of the two reactions (A+B → C, or C → A+B) takes place.
How do we decide which reaction? This is where we try to approximate reality. Here's the reasoning:
- The likelihood of the first reaction is proportional to the number of A molecules and the number of B molecules.
- This should make intuitive sense: if you have lots of both, there'll be plenty of such A's and B's combining.
- Similarly, the more we have of C, the more likely we'll have the second reaction.
Next, let
N_A = number of A molecules
N_B = number of B molecules
N_C = number of C molecules
And let
E_A = concentration of A
E_B = concentration of B
E_C = concentration of C
To get the concentration (mass per unit volume), we'll approximate the actual physical volume as:
V_total = N_A + N_B + 2*N_C
(Note: a more general and accurate formulation of volume would use the sizes of the different molecules and weight accordingly. We are simply assuming A and B molecules are similar, and C is twice as big).
Thus,
E_A = N_A / V_total
E_B = N_B / V_total
E_C = N_C / V_total
Let
R₁ = Probability that the first reaction takes place.
R₂ = Probability that the second reaction takes place.
Now, our model is going to assume that the chances of the first reaction taking place is going to be proportional to both the number of A molecules (N_A) and the number of B molecules (N_B). Thus,
R₁ is proportional to N_A N_B
Similarly,
R₂ is proportional to N_C
Because the "numbers" are proportional to concentrations, it's also true that
R₁ is proportional to E_A E_B
and
R₂ is proportional to E_C
In each case, the constant of proportionality may be slightly different. That is, it may be that the first reaction is generally more likely, all other things being equal.
To model this, we'll use constants p and q:
R₁ = E_A E_B p
R₂ = E_C q
Thus, what happens at each step of the simulation is this:
- We randomly choose between doing the first or second reaction, based on the relative proportion of R₁ and R₂.
- If it's the first, then we go ahead and remove one each of an A and a B molecule, and create a new C molecule.
- If it's the second, we remove one C molecule and create a new A and a new B molecule.

Exercise 2: Suppose you have ten A molecules and five B molecules.

Compute the concentrations (by hand) of A and B using the reasoning above.

Now run the simulation very slowly by clicking on the "Step" button. (Note: Each time you start, remember to hit the "Reset" button.)

How long before you see the second reaction taking place?
Just after you see the second reaction taking place for the first time, re-compute (by hand) the concentrations and verify that they match with the numbers in the "Stats" section.

Next, we will examine long-term behavior:

Initially, with 100 molecules of A, 50 molecules of B, and none of C, the concentrations of A, B and C are 0.667, 0.333 and 0, respectively.
After just a few steps, we can't really expect the number of A molecules to decrease drastically.
It takes a while before the system of molecules "settles down."

Exercise 3: Study the long-term behavior (starting with 100 molecules of A, 50 molecules of B):

How long (in time-steps shown) does it take for the concentrations to "settle"?
Do the concentrations ever stay fixed after a while? Why or why not?
What are the final concentrations of A, B, and C, approximately?

Let's examine how the simulation was written:

At the highest-level, the simulation looks like this:

     for n = 1 to numSteps
         draw ()                           // Draw molecules (circles)
         doReactions ()                    // Perform reaction
         collectStats ()  
     endfor
     printStats ()

The detail for each reaction:

     for n = 1 to numSteps
         // Draw molecules (circles)
         draw ()                           

         // Compute concentrations
         V = (N_A + N_B + 2*N_C)
         E_A = N_A / V
         E_B = N_B / V
         E_C = N_C / V

         // Randomly choose one reaction
         R₁ = E_A E_B p
         R₂ = E_C q
         if uniform() < R₁/(R₁ + R₂)
             // A + B → C
             Remove one A and one B, make one C
         else
             // C → A + B
             Remove one C, make one A and one C
         endif
         collectStats ()  
     endfor
     printStats ()

Exercise 4:

What is the purpose of uniform() above and how does it ensure that the reactions are selected with appropriate probabilities?
Examine the code in method standard() in MolecularSimPanel.java. What is the relevance of the spatial distance between molecules?

A more detailed model

How closely does our model correspond to reality?

In reality, molecules react when they sort-of bump into each other.
We did not consider "spatial" closeness at all.

We'll now model spatial detail:

Molecules A and B will react only if close enough.
We'll use a distance parameter that we can vary to decide what's "close enough".

Exercise 5: Use the "Model" menu and switch to the spatial model. Using the same numbers of molecules as above, see what happens when you run for about 100 time units and plot.

What do you observe? Why are there more fluctuations?
Next, look at the last row of controls. Set EndTime=100 and NumRuns=10. Then click the "simulate" button. What do you observe?
Try the same with NumRuns=100. What do you see?
For comparison, try NumRuns=10 and NumRuns=100 in the standard simulation.

On randomness in the model:

The two models so far try to simulate reality.
Sometimes, a model deliberately introduces randomness to make the model more realistic
=> e.g, the random movement of molecules.
The only source of randomness in the standard model was "which of the two reactions" to choose at each step.
What happens in the spatial model depends on the random movement of the molecules
=> there is more "chanciness" and hence more fluctuation.
Repeated trials that are then averaged can reduce the fluctations and give you the average behavior.
=> this is what we did with multiple "runs"

Here's an example to show how it works:

Concentration of A
	t=1	t=2	t=3	t=4
Run=1	0.63	0.58	0.59	0.52
Run=2	0.61	0.56	0.53	0.54
Run=3	0.62	0.59	0.55	0.51
Average	0.62	0.576	?	0.523

This shows an average of concentrations after three units of time, taken over three runs:

For each run, the concentration of A is recorded at times t=1,2,3, ....
Then, the "average run" is computed by taking the average down each column.
Thus, we get the average concentration at t=1 over all runs as 0.62.

Exercise 6: Fill in the missing number in the table above.

The spatial model is more realistic, but needs more averaging (and hence more time).
The spatial model has introduced another parameter: distance.
The spatial model does not use concentration in its model.

Exercise 7:

If the spatial model doesn't use concentration, why then is the behavior like the standard model?
Find the place in the code that has the "distance" parameter. What is the distance parameter set at? Double it and see what you get. Can you explain? What is "doubling" the distance equivalent to in terms of using other values for p and q?

A less detailed model

What we desire in a model:

Accurately matches reality.
Simple to implement, simple to understand.
Fast execution time computationally.
As little randomness as possible.
As few parameters as possible.

Consider this deterministic model:

Let
A(t) = concentration of A at time t
B(t) = concentration of B at time t
C(t) = concentration of C at time t
Now, focus on A(t):
- Consider the concentration a few moments later at time A(t+s) where s is a small value.
- In the time period, s, that has elapsed, what affects the concentration of A?
- Let's consider the change in concentration, A(t+s) - A(t).
- A(t+s) - A(t) could decrease because A molecules react with B molecules
  => the amount of this occuring depends on A(t) and B(t)
  => The more of A(t) B(t), the more that reaction will happen.
- A(t+s) - A(t) could increase because of C molecules that break down.
  => the amount of this occuring depends on C(t)
  => The more of C(t), the more this break down occurs.
- Thus, in that small time period:
  A(t+s) - A(t) decreases in proportion to A(t) B(t) (the product of these concentrations)
  A(t+s) - A(t) increases in proportion to C(t)
- Now, the total amount of change also depends on the length of the time-period, s.
  => the more s is, the more the change (because there's more time to allow more reactions).
- Putting this all together:
  A(t+s) - A(t) is proportional to s (C(t) - A(t)B(t))
We will write this as an equation:
A(t+s) - A(t) = s (C(t) - A(t)B(t))
By the same reasoning,
B(t+s) - B(t) = s (C(t) - A(t)B(t))
The change in concentration of C over this small period is different:
- C(t+s) - C(t) increases with more A(t) and more B(t).
- And it decreases in proportion to C(t).
Thus,
C(t+s) - C(t) = s (A(t)B(t) - C(t))
One minor modification:
- It may be that the the forward reaction (A+B → C) is more likely than the reverse.
- To adjust, we'll introduce multipliers for the two terms.
Thus, our revised model is:
A(t+s) - A(t) = s (K_cC(t) - K_abA(t)B(t))
B(t+s) - B(t) = s (K_cC(t) - K_abA(t)B(t))
C(t+s) - C(t) = s (K_abA(t)B(t) - K_cC(t))
This way, we can use a high value of the constant K_ab to indicate that the forward reaction is more likely, or happens at a faster rate.
Thus, we have a way to calculate the changes in concentrations as we go forward in time.
A(t+s) = A(t)+ s (K_cC(t) - K_abA(t)B(t))
B(t+s) = B(t) + s (K_cC(t) - K_abA(t)B(t))
C(t+s) = C(t) + s (K_abA(t)B(t) - K_cC(t))
How?
- Suppose we know the starting concentrations at time t=0, that is, we know A(0), B(0) and C(0).
- Suppose, further, that s=0.01, K_ab=1.0, K_c=0.5.
- We can now compute
  A(0+0.01) - A(0) = 0.01 (0.5 C(0) - 1.0 A(0)B(0))
  B(0+0.01) - B(0) = 0.01 (0.5 C(0) - 1.0 A(0)B(0))
  C(0+0.01) - C(0) = 0.01 (1.0 A(0)B(0) - 0.5 C(0))

Exercise 8: Suppose s=0.01, K_ab=1.0, K_c=0.5, and A(0)=3.0, B(0)=2.0, C(0)=1.5. Compute the concentrations of A(t), B(t), C(t) at times t=0.01, t=0.02, t=0.03 as follows:

First, calculate A(0.01), B(0.01), C(0.01).
Then, use these values to calculate A(0.02), B(0.02), C(0.02).
Then, use the values at 0.02 to calculate A(0.03), B(0.03), C(0.03).

Exercise 9: Write a Java program to compute concentrations using this model. Start by using this template. Then, modify the code to change the end-time to 1.0. What do you observe as the final concentrations of A, B and C? Plot a graph of how the concentrations change with time, and compare with previous models.

Exercise 10: How would one go about writing the code so that it's recursive? That is, our goal is to compute A(T), B(T), C(T) for some large T, which we could do recursively in terms of A(T-s), B(T-s), C(T-s), which in turn can computed from A(T-2s), B(T-2s), C(T-2s) ... etc. Write a recursive version and compare the number of arithmetic operations in the recursive approach vs. the iterative one.

How is this simulation (model) different from the others?

We have abstracted molecules away:
- We don't count numbers of molecules.
- Instead, there's this "quantity" called concentration, defined by a single number.
- This shouldn't come as a surprise because, even though we modeled individual molecules earlier, the "stats" panel and the plots showed concentrations (one number each for A, B and C).
We have abstracted away all randomness:
- No longer do we worry about how molecules move.
- The concentrations change deterministically.
This is of course a further approximation.
Are there benefits to this?
- It's faster to compute. This works for any number of molecules, including realistic values (trillions, gazillions).
- It's a simpler model. The code is only a few lines.
Any downsides?
- There are new parameters, which might be hard to measure in practice.
- There may be something to learn from the detailed molecular level that might be lost in this model.

Exercise 11: In general, how accurate should a model be? Consider this Physics experiment:

Suppose, from atop a 300-foot building on campus you fix a point on a protruding rod where you can drop a heavy (so that wind won't have an effect) object to observe where it lands on the ground.
From that same drop-point, you let down a plumb-line and mark the spot where the plumb line touches the ground.
Now you drop the heavy object.
You expect the dropped object to hit the marked spot dead on, right? Right?

Discuss with your neighboring students, whether:

The object hits the marked spot dead-on.
The object lands a little to the West.
The object lands a little to the East.

The answer to this question, it turns out, has had tremendous historical significance.

What is a differential equation?

Consider our reaction example:

Recall how we compute the concentrations evolving over time:
A(t+s) = A(t)+ s (K_cC(t) - K_abA(t)B(t))
B(t+s) = B(t) + s (K_cC(t) - K_abA(t)B(t))
C(t+s) = C(t) + s (K_abA(t)B(t) - K_cC(t))
The form of these equations is
Next value A(t+s) = current value A(t) + s * (some terms involving current values of various variables)
Equivalently, the equations specify an iterative algorithm to compute the functions A(t), B(t), C(t).
The fact that we compute in the way time evolves makes this a simulation.

Now consider the following:

Rearrange the terms in this way:
(A(t+s) - A(t)) / s = (K_cC(t) - K_abA(t)B(t))
(B(t+s) - B(t)) / s = (K_cC(t) - K_abA(t)B(t))
(C(t+s) - C(t)) / s = (K_abA(t)B(t) - K_cC(t))
This is merely a re-statement that doesn't change anything
Except: it is no longer suitable for computation.
Next, consider the limit as s → 0:
- For fixed s, (A(t+s) - A(t))/s is a function.
- Let's give these a name: Define the function A_s(t) = (A(t+s) - A(t))/s .
- Then the sequence of functions A_s possibly has a limit as s → 0.
- That limit is itself a function.
- Let's give it a name: A'(t)
  => Called the derivative (function) of function A(t)
If we do this for each of A(t), B(t), C(t), we get
A'(t) = (K_cC(t) - K_abA(t)B(t))
B'(t) = (K_cC(t) - K_abA(t)B(t))
C'(t) = (K_abA(t)B(t) - K_cC(t))
Notice the following:
- The equations specify complex relationships involving six functions
  => Three derivative functions in addition to A(t), B(t), C(t).
- These relationships are free of s.
  => We don't have to worry about whether s is small enough or too small.
- The left sides are all derivative functions.
- The right sides are all regular functions.
Such a collection of equations is called a system of differential equations.
Note: we usually call these Ordinary Differential Equations (ODEs) to differentiate them from other kinds (like partial differential equations).

What good are they?

In the days without computers, one had to reason about these systems
=> Recall, limits usually result in simplification (no s term)
ODE's are the basis for analyzing scores of real systems.
Some ODE's can be solved analytically.
Many more can be analyzed for properties (stability, long-term behavior etc) even if they can't be explicitly solved.
There are many powerful existence results that can be proved with great generality.

Relationship to computation:

Most ODE's are constructed with finite s and then taken to the limit
The equations with finite s constitute an algorithm for solving the system
=> The Euler algorithm
When the free variable t describes time, we could call this a simulation.
=> A simulation computes a system evolving over time.
Generally, for many systems, the Euler algorithm is sufficient
=> More accurate algorithms exist (e.g., the Runge-Kutta algorithm)

An unusual "reaction"

Consider the Rabbit-Lynx "reaction":

Although this is not a chemical reaction, it almost is (mathematically).
We will model a population of rabbits and lynxes living in an environment with sufficient grass for the rabbits to thrive. The lynxes, of course, eat the rabbits to thrive.
Already, we can sense the proportions involved:
- The rabbits grow at a certain rate.
- The rabbits die at a rate, proportional to how many rabbits there are and how many are being devoured by lynxes.
- The lynxes proliferate based on their consumption of rabbits, and die in numbers proportional to their population.
We will model this system:
- Let
  X(t) = the number of rabbits at time t
  Y(t) = the number of lynxes at time t
- Now consider a particular time interval from t to t+s:
- X(t+s) - X(t) is the change in the number of rabbits.
  X(t+s) - X(t) will increase from new rabbits born
  X(t+s) - X(t) will decrease from rabbits dying naturally
  X(t+s) - X(t) will decrease from rabbits killed by lynxes
  X(t+s) - X(t) is also proportional to s.
- Thus, we can write
  X(t+s) - X(t) = s * (k₁ X - k₂ X Y)
- Similarly,
  Y(t+s) - Y(t) = s * (k₂ X Y - k₃ Y)
  because the rabbits that are killed are the ones that "generate" the lynxes.
Note: we won't worry about the fact that actual numbers of rabbits are integers, and that X(t) can take a value like 0.3. Assume that we are counting in the "thousands", in which case 0.3 makes sense.
Now, we don't need anything further (other than the constants k₁, k₂, k₃) to solve this computationally.
For illustration, let us write the differential equations:
X'(t) = k₁ X - k₂ X Y
Y'(t) = k₂ X Y - k₃ Y

Exercise 12: Run the VirtualReactor program and select the "Rabbit-Lynx" model. Use these values: X(0)=1.5, Y(0)=1.0, k₁=2.4, k₂=4.2, k₃=5.1, endTime=10. (Ignore the k₄ parameter.) You should observe something quite interesting.

About the Rabbit-Lynx model:

We see that there is a natural oscillation in the populations of rabbits and lynxes:
- The populations don't ever "settle".
- Oscillation is "built in" to the model.
- We could hardly have predicted this by staring at the equations or through intuition.
The population growth of lynxes lags that of rabbits slightly as intuition would suggest.
If we think of rabbits and lynxes as interacting "molecules", the "reactions" could be written as:
grass + rabbits → 2 rabbits
rabbits + lynxes → 2 lynxes
lynxes → dead lynxes (degradation)
This model has been frequently studied, in all kinds of variations, in ecology.
The model is called a Lotka-Volterra population model.

Nonlinear chemical reactions:

Is there a chemical reaction that displays oscillation?
Indeed, there is.
The B-Z reaction :
- In 1951, the Russian chemist Belousov discovered a chemical reaction involving citric acid, acidified bromate and ceric ions displayed spectacular oscillations: turning yellow, then colorless, then yellow, in turn.
- Belousov couldn't publish his result because nobody believed it.
- Later, a student called Zhabotinsky resurrected the work and convinced others, too late unfortunately for Belousov, who died in 1970 before receiving the prestigious Lenin Prize posthumously in 1980.
See this link for some pictures of the BZ reaction.
The differential equations, after some simplification, turn out to be:
X'(t) = q Y(t) - X(t)Y(t) + X(t)(1 - X(t))
Y'(t) = (1/e) (-q Y(t) - X(t)Y(t) + f Z(t))
Z'(t) = X(t) - Z(t)
There are three variables, depicting three concentrations.

More examples

Motion along a straight line:

Let's revisit InclineSimulator.java:

	while (true) {
	    // Update time.
	    t += delT;

            // Increase velocity accordingly.
	    v = v + delT * a;

            // Increase distance according to velocity.
            d = d + delT * v;

	    if (t >= stopTime) {
                break;
	    }
	}

Note: the straight line, in this case, is along the incline.
Clearly, the differential equations are:
v'(t) = a
d'(t) = v

Exercise 13: Download and execute Missile.java. Examine the code and write down the differential equations. How many variables (functions) are involved?

Next, let's write down the differential equations for the Dubin car:

Recall: the Dubin car has two wheels, each of whose angular velocity is separately controlled.
Define the following:
(x,y) = center of the axle between the two wheels.
ω_L, ω_R = angular velocities of left and right wheels
R = radius of wheel
L = length of axle
θ = current orientation
Consider movement from center at C to C' below in a small interval of time Δ t:
- The distance CC = ½ (AA' + BB').
- Consider a small interval of time Δt. Then,
  AA' = ω_L R Δt
  BB' = ω_R R Δt
Thus,
x'(t) = R (ω_L + ω_R) cos(θ) / 2
y'(t) = R (ω_L + ω_R) sin(θ) / 2

Exercise 14: Why is this true?
The only other variable is the orientation θ(t)
=> We need a differential equation θ'(t) = ...
- Consider the case that only the right wheel moves:
- The distance moved
  AA' = ω_L R Δt
- Thus, the angle AB-AA' is:
  Δθ = ω_L R Δt / L
  
  Exercise 15: Why is this true?
- Thus,
  θ'(t) = ω_LR / L

Exercise 16: Implement a Dubin car simulator for CarSim based on these equations.

Exercise 17: What additional equation is needed for an accelerative version of the Dubin car?

Exercise 18: Write down the equations for the Simple-Car which has two controls: forward-velocity v and steering angle φ. Implement a simulator for the Simple-Car.

An example with rotational motion:

Consider a weight attached to a winch:
Define the following variables for the wheel:
R = radius of wheel
τ = torque
θ(t) = angular displacement at time t
ω(t) = angular velocity
μ(t) = angular acceleration
m_w = mass of wheel
Define the following variables for the cable:
T = tension in the cable
Define the following variables for the load:
m = mass of load
y(t) = vertical displacement of load
v(t) = vertical velocity of load
a(t) = vertical acceleration of load
Now apply Newton's laws to each component. Let's start with the wheel:
- The wheel has torque generated from the winch's motor: τ
- The tension in the cable pulls the other way, and creates a reverse torque of TR
- Thus, the difference is what determines the angular motion of the wheel:
  m_w R² μ = τ - TR
  Note: we've used moment of inertia (mass analog for angular motion)
Now the load:
- The upward force is the cable's tension: T
- The downward force is the weight mg
- The difference explains the motion:
  ma = T - mg
What are the differential equations?
- First, the obvious ones relating acceleration, velocity and displacement for the wheel:
  ω'(t) = μ(t)
  θ'(t) = ω(t)
- The same for the load:
  v'(t) = a(t)
  y'(t) = v(t)
Observe that a(t) = Rμ(t)

Exercise 19: Why?
Finally, let's relate these through the tension:
- From the wheel:
  T = (τ - m_wR²μ) / R
- From the load:
  T = mRμ - mg
- Equating the two and solving for μ:
  μ = (τ - mgR) / (mR²+m_wR²)
Thus, the differential equations that describe the system are:
μ(t) = (τ - mgR) / (mR²+m_wR²)
ω'(t) = μ(t)
θ'(t) = ω(t)
y(t) = Rθ(t)

Exercise 20: Download and examine Winch.java.

Identify the code corresponding to the differential equations we derived earlier.
Handcode different values of torque to see what value is sufficient to lift the load. Is there a value that will lift the load and yet avert a collision with the winch?

Exercise 21: Set isVertical=false in the program and execute. You should see a version without gravity (as if the load were on a frictionless surface). What value of torque is sufficient to pull the load?