In the process of creating a physics simulation we start by inventing a mathematical model and finding the differential equations that embody the physics. The next step is getting the computer to solve the equations, a process that goes by the name numerical analysis.
For simple models you can use calculus, trigonometry, and other math techniques to find a function which is the exact solution of the differential equation. This is called the analytic solution (because you use analysis to figure it out). It is also referred to as a closed form solution. (BTW, college classes on differential equations are all about finding analytic solutions).
An analytic solution is preferred because you can get some insight from the equation. For example, in the analytic solution of the Single Spring simulation we can see how mass m or spring stiffness k affects the frequency of oscillation, just by looking at the equation:
$$x(t) = x_0 \cos(\sqrt{k/m} \; t)$$
In contrast, with a numerical solution you would have to run lots of experiments to figure out relationships like that; just like in real life you would set up experiments with different masses or springs, make measurements, and then make a graph that shows how frequency varies for different masses or spring stiffness.
Most physics simulations are too complicated to be able to find an analytic solution. Instead we use Numerical Analysis to find an approximate numerical solution, which is just a list of numbers. The numbers are model's variables at each moment in time. For example, here are the first few moments of the Single Spring simulation
t x v
0.00000 −2.00000 0.00000
0.02500 −1.99626 0.29906
0.05000 −1.98507 0.59552
0.07500 −1.96651 0.88827
0.10000 −1.94070 1.17623
0.12500 −1.90775 1.45837
0.15000 −1.86783 1.73364
0.17500 −1.82113 2.00105
0.20000 −1.76785 2.25965
0.22500 −1.70823 2.50851
0.25000 −1.64252 2.74675
This list of numbers is what drives the visual representation of the simulation: for the current value of time t we move the block to the position given by the x variable at that time.
How does it work? How do we generate that list of numbers?
We start with the initial values of the variables, which is the first row in the table above. We then use the differential equations for the simulation to calculate the variables after a very brief period of time. The differential equation tells us the rate of change in each variable at any point in time. By adding together small changes over time (also known as integrating over time) we can predict the future.
For the Single Spring simulation the differential equations are:
| x' = v | (1) |
| v' = −k⁄m x − b⁄m v | (2) |
The parameter values for spring stiffness k , mass m , and damping b are:
k = 3.0
m = 0.5
b = 0.1
With these parameters, equations (1) and (2) simplify to
| x' = v | (3) |
| v' = −6 x − 0.2 v | (4) |
The initial conditions at time t = 0 are
x = −2.0
v = 0.0
We can find the rate of change of the variables at time t = 0 by substituting the above values into equations (3) and (4)
x' = 0
v' = −6 (−2) − 0.2 (0) = +12
Since we know the starting values of the variables and their rates of change, we can predict the future values. At least for a very short amount of time in the future, assuming that the rates of change do not themselves change too much!
To show how numerical integration works, we use the simplest possible method: Euler's method. It is simple but not very accurate. We'll mention more accurate methods below.
In solving differential equations (3) and (4), assume we know the state of the system at time tn to have position xn and velocity vn .
To get the state of the system at time tn+1 = tn + Δt after a small amount of time Δt has elapsed:
| xn+1 = xn + Δt vn | (5) |
| vn+1 = vn + Δt (− 6 xn − 0.2 vn) | (6) |
Equations (5) and (6) use the rates given by the differential equations (3) and (4) to make a crude estimate of the new value after the small amount of time Δt has elapsed.
Now we can use equations (5) and (6) to make the numerical predictions. We start at time t = 0 with
x0 = −2
v0 = 0
With n = 0 and a small timestep of Δt = 0.025 , the equations (5) and (6) give the state at time t1 = t0 +Δt = 0.025
x1 = x0 + 0.025 v0 = −2 + 0.025 * 0 = −2
v1 = v0 + 0.025 (−6 x0 − 0.2 v0) = 0 + 0.025 (−6 (−2) − 0.2 (0)) = 0.3
With n = 1 in equations (5) and (6) we can get the state at time t2 = t1 + Δt = 0.050
x2 = x1 + 0.025 v1 = −2 + 0.025 (0.3) = −1.9925
v2 = v1 + 0.025 (−6 x1 − 0.2 v1) = 0.3 + 0.025 (−6 (−2) − 0.2 (0.3)) = 0.5985
With n = 2 in equations (5) and (6) we can get the state at time t3 = t2 + Δt = 0.075
x3 = x2 + 0.025 v2 = −1.9925 + 0.025 (0.5985) = −1.9775375
v3 = v2 + 0.025 (−6 x2 − 0.2 v2)
= 0.5985 + 0.025 (−6 (−1.9925) − 0.2 (0.5985)) = 0.8943825
Our approximate numerical solution thus far:
t x v
0.00000 −2.00000 0.00000
0.02500 −2.00000 0.30000
0.05000 −1.99250 0.59850
0.07500 −1.97754 0.89438
Here again for comparison is the solution that is calculated in the Single Spring simulation using the highly accurate Runge-Kutta Method
t x v
0.00000 −2.00000 0.00000
0.02500 −1.99626 0.29906
0.05000 −1.98507 0.59552
0.07500 −1.96651 0.88827
Our simple-minded solution isn't that far off! But if we compare this to a real world example of an actual spring and mass, we would see that this method is very inaccurate. Or if we look at whether energy is being conserved we will also see it has poor accuracy. For many simulations, this method is so bad that it becomes unstable and will "blow up" with solutions going to infinity.
There are several methods for integrating a differential equation over time. They have various tradeoffs of accuracy vs. computational cost. The Wikipedia article Numerical methods for ordinary differential equations describes several of them.
The methods used in the myphysicslab simulations are:
These numerical methods are defined with equations that can get a bit confusing. To get used to the notation, here is how the simple Euler method we used above could be notated. In more general form we can write equations (5) and (6) like this:
| xn+1 = xn + Δt f (xn, vn) | (7) |
| vn+1 = vn + Δt g(xn, vn) | (8) |
where the functions f, g are from the differential equations that we are solving written like this:
| x' = f (x, v) | (9) |
| v' = g(x, v) | (10) |
The problem with the simple Euler method numerical solver is that it only uses the rate of change at the beginning of the time period from tn to tn + Δt . Yet that rate of change (the derivative) is itself changing over this period! The more accurate numerical methods estimate how the derivatives in equations (9) and (10) change over the time period.
