CS 4: Lecture 7
Wednesday, February 8, 2006
PHYSICAL SIMULATIONS
====================
In lab, we've seen how to compute trajectories in which a projectile undergoes
a constant acceleration, using the equation x = x0 + v0 t + a t^2 / 2. But
what if the acceleration is not constant? The equation doesn't work any more.
If we know the velocity or acceleration as a function of time, there are two
ways we might be able to determine the position of an object at a given time.
- An _analytical_solution_: find an explicit expression for the position as a
function of time, like x(t) = x0 + v0 t + a t^2 / 2. Unfortunately, this is
not always possible.
- A _numerical_solution_: find an _approximate_ solution based on a computer
simulation. Unfortunately, this doesn't give 100% accuracy, but we can often
get very close. This is what we'll look at today.
Creating a computer simulation takes several steps.
- Physics: understand how the quantities are related.
- Discretization: the position of an object is a continuous function of time,
but computers can only do a finite number of computations, so we're going to
look at a discrete (finite) set of "points in time." We need to figure out
an _algorithm_ to do this.
- Simulation: write and run code that calculates the discretized values.
Understand the output.
Let's look at these steps.
Physics
-------
Suppose you drive a car whose velocity at time t is v(t) = sin t miles/minute,
where t is measured in minutes, and the sine is taken over degrees (not
radians). You start at time t = 0 with a velocity of v(0) = 0, and stop when
your velocity peaks at 1 mile/minute after 90 minutes. How far have you
driven?
Let x(t) express how far you've driven at time t. Remember that velocity is
the derivative of distance with respect to time, so
dx(t) o <-- A reminder that t is in degrees, not radians.
----- = v(t) = sin t .
dt
This is called a _differential_equation_. To solve it, you also need to know
the _initial_condition_ x(0) = 0. (If you started driving at x(0) = 10, you'd
obviously end up in a different place than if you start driving at x(0) = 0.)
The analytical way to solve this differential equation is to integrate, giving
/\ t pi o
x(t) = | sin --- y dy. <-- Here we convert y to radians so you
\/ 0 180 can integrate the sine function.
However, there are many differential equations that arise in practice that
nobody knows how to integrate or solve. For example, if you want to account
for the effects of wind resistance when you simulate slingshotting a ball along
a trajectory, there is no known analytical solution. If you want to simulate
three planets orbiting each other under gravitational forces, there is no known
analytical solution. (There's a solution for two planets--they revolve each
other in elliptical orbits--but not for three.) So let's look at how to
approximate the integral numerically.
Discretization
--------------
Here's the idea. Suppose we know the position x(t) at time t, and we want to
know the position x(t + T) at time t + T, where T is a small positive number.
Suppose that the velocity is a constant, v, during this period of time. Then
x(t + T) = x(t) + v T.
Unfortunately, in our problem, v(t) is not constant! But if we choose T to be
really small, then v(t) is _close_ to constant during that period of time.
So although we won't get x(t + T) exactly right, we'll be close. And the
smaller we make T, the less v(t) will vary during that time, so the closer
we'll be. If we could make T as small as the infinitesimals used in calculus,
we'd get exactly the right answer. Computers can't do infinitesimals, but they
can make T pretty small if you're willing to wait long enough.
Here's a mathematical justification for this approximation. If we express
x(t + T) as a Taylor series, expanded around t, we have
2 3
dx(t) 1 d x(t) 2 1 d x(t) 3
x(t + T) = x(t) + ----- T + - ------ T + - ------ T + ...
dt 2 2 6 3
dt dt
If T is small enough, then T^2 is really small, and T^3 is absolutely tiny.
The smaller you make T, the more the first two terms dominate all the rest.
So if T is small enough, we can ignore the terms involving T^2, T^3, etc.,
giving the approximation
x(t + T) ~ x(t) + v(t) T.
Our strategy is to approximate a sequence of values of x(t), each at a
different time t.
x(0) = 0.
x(T) ~ x(0) + v(0) T. <- The "~" means "is approximately equal"
x(2T) ~ x(T) + v(T) T. and is actually written as two squiggles,
x(3T) ~ x(2T) + v(2T) T. one over the other, when I'm not stuck
x(4T) ~ x(3T) + v(3T) T. with ASCII.
... and so on, until we get to the end of the 90-minute trip. To compute the
distance x(90), repeat 90/T times.
This is called a _discretization_ of the original integral we wanted to solve,
because we compute a finite number of discrete _timesteps_ of duration T,
instead of an infinite number of infinitesimal timesteps. This particular
discretization is called _Euler's_method_, or alternatively
_forward_differencing_. We'll see a more accurate discretization soon.
This strategy is also called _numerical_integration_, because we are
approximating the integral numerically. We are approximating the area under
the sine curve by adding the areas of a bunch of tall, thin rectangles (one for
each timestep).
____
__||||
__||||||
_|||||||| Each rectangle has width T and height v(t), where t is
_||||||||| the time at the left edge of the rectangle.
_||||||||||
Of course, the flat tops of the rectangles don't match the curve of the sine,
so our approximation of x(t) is not perfectly accurate. But we can improve the
accuracy by making T smaller. However, we have to do 90/T timesteps, so better
accuracy means a longer computation time.
Simulation
----------
Now we can write code to compute the length of our drive.
public static void main(String[] args) {
double duration = 90.0; // Duration of simulation
long timesteps = 10; // Number of timesteps
double steptime = duration / (double) timesteps; // Time for one timestep
double distance = 0.0; // Distance traveled
long step = 0; // Current step number
while (step < timesteps) {
double time = steptime * (double) step;
/* We must convert `time' from degrees to radians to use Math.sin */
double velocity = Math.sin(Math.PI / 180.0 * time);
distance = distance + velocity * steptime;
step++;
}
System.out.println("You drove " + distance + " miles.");
}
Each iteration of the loop performs one timestep of the simulation.
Let's look at the output of this program...for a bunch of different choice of
the variable `timesteps'.
# of timesteps | Miles driven (output)
-----------------|------------------------
1 | 0.
2 | 31.
5 | 47.8
10 | 52.6
100 | 56.84
1,000 | 57.251
10,000 | 57.2912
100,000 | 57.29533
As the number of timesteps increases, the solutions _converge_--you can see
they're getting closer to the same answer. In this case, we can evaluate the
integral exactly, so we know the exact answer is about 57.295779. Not only is
the simulation converging--it's converging to the right answer.
(Unfortunately, there are many problems in scientific computing where it's
quite hard to devise a numerical method that converges to the right answer.)
A Better Discretization
-----------------------
Instead of using Euler's method, let's use a method called the _trapezoid_rule_
to integrate v(t). Instead of approximating the integral with a bunch of
rectangles, we'll use a bunch of tall thin trapezoids whose tops are tilted.
Both endpoints of each trapezoid top lie on the sine curve.
_-/
_-/|||
_/||||||
/||||||||
/|||||||||
/||||||||||
Mathematically, this means that during the timestep that takes us from t to
t + T, we don't assume that the velocity is v(t) over the whole interval.
Instead, we assume the velocity is the _average_ of v(t) and v(t + T).
v(t) + v(t + T)
x(t + T) ~ x(t) + --------------- T.
2
We can implement this method by changing just one line of code.
double velocity = 0.5 * Math.sin(Math.PI / 180.0 * time) +
0.5 * Math.sin(Math.PI / 180.0 * (time + steptime));
Let's look at the output using our improved program.
# of timesteps | Euler's method | Trapezoid rule
-----------------|--------------------------------------
1 | 0. | 45.
2 | 31. | 54.3
5 | 47.8 | 56.82
10 | 52.6 | 57.18
100 | 56.84 | 57.2946
1,000 | 57.251 | 57.295768
10,000 | 57.2912 | 57.29577939
100,000 | 57.29533 | 57.2957795119
The exact answer is about 57.2957795131. Observe that the trapezoid rule
converges to the correct solution a _lot_ faster than Euler's method! In 1,000
timesteps, the trapezoid rule is already more accurate than Euler's method
after 100,000 timesteps.
Postscript: Why the Trapezoid Rule Works (not examinable)
----------------------------------------------------------
Why is the trapezoid rule so much more accurate? Expand x'(t + T) around t.
2 3
dx(t + T) dx(t) d x(t) 1 d x(t) 2
--------- = ----- + ------ T + - ------ T + ...
dt dt 2 2 3
dt dt
Take the Taylor series for x(t + T), add T/2 times the left-hand side of this
equation, and subtract T/2 times the right-hand side.
3
1 dx(t) 1 dx(t + T) 1 d x(t) 3
x(t + T) = x(t) + - ----- T + - --------- T - -- ------ T + ...
2 dt 2 dt 12 3
dt
Ta-daa! The T^2 terms magically cancelled each other out. The trapezoid
method uses the first three terms, so its error is proportional to T^3. If T
is small, then T^3 is tiny--smaller than the T^2 error of Euler's method.