Simulation of the physical world

How would you approach the simulation of rain, well, or any other lengthy physical process?

Simulation, be it rain, air flow over the wing of the plane or falling down the steps of the link (remember the rainbow spring toy from childhood?), You can imagine if you know the following:
')

1. The state of everything at the time of the start of the simulation.
2. How does this state change from one time point to another?

By “state of everything,” I mean any variable data that either determines how the environment looks, or changes that occur over time. As an example of a state, one can cite the position of a raindrop, the direction of the wind, or the speed of each part of the spring of a link.

If we assume that all our state of the environment is one big vector $$$\ vec {y}$ , then we can formulate the data we need, mentioned above, in the following:

1. Find value $$$y_0$ satisfying $$$y (t_0) = y_0$ ?
2. Find function $$$f$ such that $$$\ frac {dy (t)} {dt} = f (t, y (t))$ .

Why do we need to store the state of all in one vector, I will explain a little later. This is one of those cases where it seems that we are overdoing the task summary, but I promise there are some interesting things about this approach.

Now let's see how you can store all the information about the environment in one vector on a simple example. Suppose we have 2 objects in a 2D simulation. Each object is determined by its position. $$$\ vec x$ and speed $$$\ vec v$ .


So, to get a vector $$$\ vec y$ we need to combine vector $$$\ vec x_1, \ vec v_1, \ vec x_2 ~ and ~ \ vec y_2$ into one, consisting of 8 components, like this:

$$

$$\ vec y = \ begin {bmatrix} \ vec x_1 \\ \ vec v_1 \\ \ vec x_2 \\ \ vec v_2 \ end {bmatrix} = \ begin {bmatrix} x_ {1x} \\ x_ {1y} \ \ v_ {1x} \\ v_ {1y} \\ x_ {2x} \\ x_ {2y} \\ v_ {2x} \\ v_ {2y} \ end {bmatrix}$$

If you are confused, why do we want to find $$$f (t, y (t))$ , not the initial definition $$$\ frac {dy (t)} {dt}$ , the meaning is that we need the derivative as an expression depending only on the current state, $$$y (t)$ , constants and time itself. If this is not possible, then surely any information about the state was not taken into account.

Initial state

To determine the initial state of the simulation, you need to define a vector $$$\ vec y_0$ . Thus, if the initial state of our example with two objects looks like this:


Then in a vector form it can be represented as follows:

$$

$$\ vec x_1 (t_0) = \ begin {bmatrix} 2 \\ 3 \ end {bmatrix}, \ vec v_1 (t_0) = \ begin {bmatrix} 1 \\ 0 \ end {bmatrix}, \ vec x_2 (t_0 ) = \ begin {bmatrix} 4 \\ 1 \ end {bmatrix}, \ vec v_2 (t_0) = \ begin {bmatrix} -1 \\ 0 \ end {bmatrix}$$

If we combine it all into one vector, we get the one we need. $$$\ vec y_0$ :

$$

$$\ vec y_0 = \ vec y (t_0) = \ begin {bmatrix} x_ {1x} \\ x_ {1y} \\ v_ {1x} \\ v_ {1y} \\ x_ {2x} \\ x_ {2y } \\ v_ {2x} \\ v_ {2y} \ end {bmatrix} = \ begin {bmatrix} 2 \\ 3 \\ 1 \\ 0 \\ 4 \\ 1 \\ -1 \\ 0 \ end { bmatrix}$$

Derivative function

$$$\ vec y_0$ determines the initial state and now all we need is to find a way to transition from the initial state to the state what is happening moment after it, and from the obtained state a little further in time and so on.

To do this, we solve the equation $$$\ frac {dy (t)} {dt} = f (t, y (t))$ for $$$f$ . Find the derivative of $$$y (t)$ :

$$

$$\ frac {dy (t)} {dt} = \ frac {d} {dt} \ begin {bmatrix} x_ {1x} \\ x_ {1y} \\ v_ {1x} \\ v_ {1y} \\ x_ {2x} \\ x_ {2y} \\ v_ {2x} \\ v_ {2y} \ end {bmatrix} = \ begin {bmatrix} \ frac {dx_ {1x}} {dt} \\ \\ \ frac {dx_ {1y}} {dt} \\ \\ \ frac {dv_ {1x}} {dt} \\ \\ \ frac {dv_ {1y}} {dt} \\ \\ \ frac {dx_ {2x} } {dt} \\ \\ \ frac {dx_ {2y}} {dt} \\ \\ \ frac {dv_ {2x}} {dt} \\ \\ \ frac {dv_ {2y}} {dt} \ end {bmatrix}$$

Wow! High turned formula! But you can bring it into a more readable form, if we divide our vector $$$\ vec y$ back to the component parts.

$$

$$\ begin {aligned} \ frac {d \ vec x_1 (t)} {dt} = \ begin {bmatrix} \ frac {dx_ {1x}} {dt} \\ \\ \ frac {dx_ {1y}} { dt} \ end {bmatrix}, \ frac {d \ vec v_1 (t)} {dt} = \ begin {bmatrix} \ frac {dv_ {1x}} {dt} \\ \ frac {dv_ {1y} } {dt} \ end {bmatrix} \\ \\ \ frac {d \ vec x_2 (t)} {dt} = \ begin {bmatrix} \ frac {dx_ {2x}} {dt} \\ \\ \ frac {dx_ {2y}} {dt} \ end {bmatrix}, \ frac {d \ vec v_2 (t)} {dt} = \ begin {bmatrix} \ frac {dv_ {2x}} {dt} \\ \\ \ frac {dv_ {2y}} {dt} \ end {bmatrix} \ end {aligned}$$

$$$\ vec x_1$ and $$$\ vec x_2$ bound by similar rules, as well as $$$\ vec v_1$ with $$$\ vec v_2$ , so despite a bunch of expressions on top, all we really need to find are the following 2 things:

$$

$$\ frac {d \ vec x} {dt} ~~ \ text {and} ~ \ frac {d \ vec v} {dt}$$

The quality of the simulation depends on the definition of these two derivatives, it is in them that the whole force lies. And so that the simulation does not resemble a program where everything happens in a chaotic way, you can turn to physics for inspiration.

Kinematics and Dynamics

Kinematics and dynamics - the necessary ingredients to create an interesting simulation. Let's start with the basics in our example.

Responsible for the position in space $$$\ vec x$ , and the first derivative of the position of a point in time is its speed $$$\ vec v$ . In turn, the derivative of speed over time is acceleration, $$$\ vec a$ .

It may seem that we have already found the function we need. $$$f$ because we already know the following:

$$

$$\ frac {d \ vec x} {dt} = \ vec v ~~ \ text {and} ~ \ frac {d \ vec v} {dt} = \ vec a$$

And in fact, we coped brilliantly with $$$\ frac {d \ vec x} {dt}$ because $$$\ vec v$ this is part of our state vector $$$\ vec y (t)$ but we need a little more to disassemble the second formula, because with $$$\ vec a$ not so simple.

Here Newton's second law will help us: $$$\ vec F = m \ vec a$ . If we assume that the mass of our objects is known, then by rearranging the variables in the expression, we get:

$$

$$\ frac {d \ vec v} {dt} = \ vec a = \ frac {\ vec F} {m}$$

So, wait a minute, $$$\ vec a$ and $$$\ vec F$ not part $$$\ vec y (t)$ therefore it is difficult to call it promotion (remember, we need a derivative function depending only on $$$\ vec y (t)$ and $$$t$ ). Nevertheless, we have moved forward because we have found all the useful formulas that determine the behavior of objects in our physical world.

Suppose that in our simple example, the only force that affects objects is gravitational attraction. In that case, we can determine $$$\ vec F$ using the law of world wide Newton:

$$

$$F = G \ frac {m_1 m_2} {r ^ 2}$$

Where $$$G$ it's a gravitational constant $$$6.67 \ times 10 ^ {- 11} \ frac {Nm ^ 2} {kg ^ 2}$ , but $$$m_1$ and $$$m_2$ these are the masses of our objects (which, we assume, are constants).

To create the simulation itself, we also need the direction and how to specify $$$r$ through vector components $$$\ vec y (t)$ . Assuming that $$$\ vec F_1$ this is the force acting on the first object, then you can do it as follows:

$$

$$\ begin {aligned} \ vec F_1 & = G \ frac {m_1 m_2} {| \ vec x_2 - \ vec x_1 | ^ 2} \ left [\ frac {\ vec x_2 - \ vec x_1} {| \ vec x_2 - \ vec x_1 |} \ right] = G \ frac {m_1 m_2 (\ vec x_2 - \ vec x_1)} {| \ vec x_2 - \ vec x_1 | ^ 3} \\ \\ \ vec F_2 & = G \ frac {m_2 m_1} {| \ vec x_1 - \ vec x_2 | ^ 2} \ left [\ frac {\ vec x_1 - \ vec x_2} {| \ vec x_1 - \ vec x_2 |} \ right] = G \ frac {m_2 m_1 (\ vec x_1 - \ vec x_2)} {| \ vec x_1 - \ vec x_2 | ^ 3} \ end {aligned}$$

To summarize The change of states in our system of two objects is fully expressed in terms of variables. $$$\ vec x_1, \ vec v_1, \ vec x_2, ~ and ~ \ vec v_2$ . And the changes are:

$$

$$\ begin {aligned} \ frac {d \ vec x_1} {dt} & = \ vec v_1 \\ \\ \ frac {d \ vec x_2} {dt} & = \ vec v_2 \\ \\ \ frac {d \ vec v_1} {dt} & = \ vec a_1 = \ frac {\ vec F_1} {m_1} = G \ frac {m_2 (\ vec x_2 - \ vec x_1)} {| \ vec x_2 - \ vec x_1 | ^ 3} \\ \\ \ frac {d \ vec v_2} {dt} & = \ vec a_1 = \ frac {\ vec F_1} {m_1} = G \ frac {m_1 (\ vec x_1 - \ vec x_2)} { | \ vec x_1 - \ vec x_2 | ^ 3} \ end {aligned}$$

Now we have everything that distinguishes our simulation from all other simulations: $$$\ vec y_0$ and $$$f$ .

$$

$$\ begin {aligned} \ vec y_0 & = \ vec y (0) & = \ begin {bmatrix} \ vec x_1 (0) \\ \ vec v_1 (0) \\ \ vec x_2 (0) \\ \ vec v_2 (0) \ end {bmatrix} & = \ begin {bmatrix} (2, 3) \\ (1, 0) \\ (4, 1) \\ (-1, 0) \ end {bmatrix} \\ \\ f (t, y (t)) & = \ frac {d \ vec y} {dt} (t) & = \ begin {bmatrix} \ frac {d \ vec x_1} {dt} (t) \\ \\ \ frac {d \ vec v_1} {dt} (t) \\ \\ \ frac {d \ vec x_2} {dt} (t) \\ \\ \ frac {d \ vec v_2} {dt} ( t) \ end {bmatrix} & = \ begin {bmatrix} \ vec v_1 (t) \\\\ G \ frac {m_2 \ big (\ vec x_2 (t) - \ vec x_1 (t) \ big)}} | \ vec x_2 (t) - \ vec x_1 (t) | ^ 2} \\ \\ \ vec v_2 (t) \\ \\ G \ frac {m_1 \ big (\ vec x_1 (t) - \ vec x_2 (t) \ big)} {| \ vec x_1 (t) - \ vec x_2 (t) | ^ 2} \ end {bmatrix} \ end {aligned}$$

But how, having a strictly defined simulation, turn it into a beautiful animation?

If you have had experience writing a simulation or game, then maybe you offer something like this:

x += v * delta_t v += F/m * delta_t 

But let's stop a little and see why it works.

Differential equations

Before we begin to implement, you need to decide what information we already have and what we need. We have value $$$y_0$ which satisfies $$$y (t_0) = y_0$ there is also $$$f$ satisfying $$$\ frac {dy} {dt} (t) = f (t, y (t))$ . And we need a function that will give us the state of the system at any time. Mathematically, we need a function. $$$y (t)$ .

Having this in mind and looking closer $$$\ frac {dy} {dt} (t) = f (t, y (t))$ You can notice that this equation links $$$y$ with its derivative $$$\ frac {dy} {dt}$ . This means that our equation is differential! Ordinary first-order differential equation, to be more precise. If solved, we will find the function $$$y (t)$ .

Task of finding $$$y (t)$ according to $$$y_0$ and $$$f$ called the Cauchy problem .

Numerical integration

For some examples of Cauchy problems, you can easily find the answer by an analytical method, but in complex simulations an analytical approach can be very difficult. Therefore, we will try to find a way to find an approximated solution of the problem.

For example, take the simple Cauchy problem.
Given: $$$y (0) = 1$ and $$$\ frac {dy} {dt} (t) = f (t, y (t)) = y (t)$ . Find an approximate solution for $$$y (t)$ .

Consider the problem from a geometric point of view and look at the value and tangent at the point $$$t = 0$ . From what is given to us, we have $$$y (0) = 1$ and $$$\ frac {dy} {dt} (t) = y (t) = 1$


We don't know what it looks like yet. $$$y (t)$ but we know that near the point $$$t = 0$ , value close to tangent. Now we will try to calculate $$$y (0 + h)$ for small value $$$h$ taking advantage of the tangent. First try $$$h = 0.5$ .


If we paint, then we approximate the value $$$y (h)$ in the following way:

$$

$$\ begin {aligned} y (h) \ approx y (0) + h \ frac {dy} {dt} (0) & = y (0) + hf (0, y (0)) \\ & = y (0) + hy (0) \\ & = 1 + h \ end {aligned}$$

So for $$$h = 0.5, y (h) \ approx 1.5$ .
Go
Now we can continue to calculate for other points. Although, of course, we found not the exact value $$$y (h)$ , but if our approximate value is very close to exact, then the approximated tangent will also be very close to real!

$$$$ display $$ \ begin {aligned} f (t, y (t)) & = y (t) \\ f (0.5,1.5) & = 1.5 \ end {aligned} $$ display $$

Next, let's move on to $$$h$ units to the right tangent.


Repeat the process and get the angular coefficient of the tangent $$$f (t, y (t)) = f (1,2.25) = 2.25$ :


The procedure can be performed recursively and for this we derive the formula:

$$

$$\ begin {aligned} t_ {i + 1} & = t_i + h \\ y_ {i + 1} & = y_i + h f (t_i, y_i) \ end {aligned}$$

This numerical method for solving differential equations is called the Euler method. For the general case, the step is x += v * delta_t .

In our particular case, the step by step solution looks like this:

$$

$$\ begin {aligned} t_ {i + 1} & = t_i + h \\ y_ {i + 1} & = y_i + h y_i \ end {aligned}$$

Using this method, it is convenient to present the results in the form of a table:

$$

$$\ begin {aligned} t_0 & = 0, & y_0 & = 1 & & & \\ t_1 & = 0.5, & y_1 & = y_0 + h y_0 & = & 1 + 0.5 (1) & = & 1.5 \\ t_2 & = 1, & y_2 & = y_1 + h y_1 & = & 1.5 + 0.5 (1.5) & = & 2.25 \\ t_3 & = 1.5, & y_3 & = y_2 + h y_2 & = & 2.25 + 0.5 (2.25) & = & 3.375 \\ \ end {aligned}$$

It turns out that our problem has a beautiful analytical solution. $$$y (t) = e ^ {t}$ :


What do you think will happen if you reduce the step in the Euler method?


The difference between the approximated and exact solutions decreases with decreasing $$$h$ ! In addition, in addition to the step reduction, you can use other methods of numerical integration, which can lead to a better result, such as the method of averages rectangles , the Runge-Kutta method and the Adams method .

It's time to code!

With the same success as we derived a mathematical representation of a description of a simulation, we can write the implementation of a simulation programmatically.

Since I'm most familiar with JavaScript, and I like the clarity that annotations add to the code; all the examples will be written in TypeScript .

And we will start with the version in which it meant that $$$y$ This is a one-dimensional array of numbers, just like in our mathematical model.

 function runSimulation( // y(0) = y0 y0: number[], // dy/dt(t) = f(t, y(t)) f: (t: number, y: number[]) => number[], //     render: (y: number[]) => void ) { //    1/60    //    60fps         const h = 1 / 60.0; function simulationStep(ti: number, yi: T) { render(yi) requestAnimationFrame(function() { const fi = f(ti, yi) // t_{i+1} = t_i + h const tNext = ti + h // y_{i+1} = y_i + hf(t_i, y_i) const yNext = [] for (let j = 0; j < y.length; j++) { yNext.push(yi[j] + h * fi[j]); } simulationStep(tNext, yNext) } } simulationStep(0, y0) } 

It is not always convenient to operate with one-dimensional arrays; you can abstract the functions of addition and multiplication of the simulation process into an interface and get a brief generalized implementation of the simulation using TypeScript Generics .

 interface Numeric<T> { plus(other: T): T times(scalar: number): T } function runSimulation<T extends Numeric<T>>( y0: T, f: (t: number, y: T) => T, render: (y: T) => void ) { const h = 1 / 60.0; function simulationStep(ti: number, yi: T) { render(yi) requestAnimationFrame(function() { // t_{i+1} = t_i + h const tNext = ti + h // y_{i+1} = y_i + hf(t_i, y_i) const yNext = yi.plus(f(ti, yi).times(h)) simulationStep(yNext, tNext) }) } simulationStep(y0, 0.0) } 

The positive side of this approach is the ability to concentrate on the basis of a simulation: what exactly this simulation distinguishes from any other. We use the simulation example with the two objects mentioned above:


If you podshamanit with numbers, then you can get a simulation of the moon's orbit!

Simulation of the orbit of the moon, 1 pix. = 2500 km. 1 sec Simulation is 1 day on Earth. The proportion of the moon to the earth is increased 10 times

Collisions and limitations

The above mathematical model really simulates the physical world, but in some cases the numerical integration method, unfortunately, breaks down.

Imagine a simulation of a ball jumping on the surface.

The simulation state can be described as:

$$

$$\ vec y = \ begin {bmatrix} x \\ v \ end {bmatrix}$$

Where $$$x$ this is the height of the ball above the surface, and $$$v$ his speed If you release the ball from a height of 0.8 meters, you get:

$$

$$\ vec y_0 = \ begin {bmatrix} 0.8 \\ 0 \ end {bmatrix}$$

If you draw a graph $$$x (t)$ , then we get something like this:


During the fall of the ball, the derivative function $$$f$ it is calculated quite easily:

$$

$$f (t, y (t)) = \ frac {dy} {dt} (t) = \ begin {bmatrix} \ frac {dx} {dt} (t) \\ \\ \ frac {dv} {dt } (t) \ end {bmatrix} = \ begin {bmatrix} v \ \ a \ end {bmatrix}$$

With the acceleration of free fall, $$$a = -g = -9.8 \ frac {m} {s ^ 2}$ .

But what happens when the ball touches the surface? The fact that the ball has reached the surface can be recognized by $$$x = 0$ . But with numerical integration, at one point in time the ball can be above the surface, and already at the next moment below it: $$$x_i> 0, x_ {i + 1} <0$ .

It would be possible to solve this problem by determining the moment of collision. $$$t_ {c} ~ (t_i <t_c <t_ {i + 1})$ . But even if this moment is found, how to determine the acceleration $$$\ frac {dv} {dt}$ so that it changes in the opposite direction.

You can, of course, define a collision in a limited period of time and apply another force $$$F$ for this length of time $$$\ Delta t$ , but it is much easier to define a discrete constant bounding simulation.

And in order to reduce the amount of surface penetration by the ball, it is possible to calculate several simulation steps at once for one tick. Together with this, our simulation code will change like this:

 function runSimulation<T extends Numeric<T>>( y0: T, f: (t: number, y: T) => T, applyConstraints: (y: T) => T, iterationsPerFrame: number, render: (y: T) => void ) { const frameTime = 1 / 60.0 const h = frameTime / iterationsPerFrame function simulationStep(yi: T, ti: number) { render(yi) requestAnimationFrame(function () { for (let i = 0; i < iterationsPerFrame; i++) { yi = yi.plus(f(ti, yi).times(h)) yi = applyConstraints(yi) ti = ti + h } simulationStep(yi, ti) }) } simulationStep(y0, 0.0) } 

And now you can write the code of our bouncing ball:

Attention developers!

Although this model has its advantages, it does not always lead to productive simulations. For me, such a framework is useful for representing the behavior of a simulation, even if a lot of extra things happen in it.

For example, a rain simulation at the beginning of this article [ approx. In the original article, at the beginning a beautiful interactive simulation of rain was inserted, I recommend to see firsthand] it was not created using the template described in the article. It was an experiment using Entity – component – ​​system . Simulation sources can be found here: rain simulation on GitHub .

See you later!

I find the intersection of mathematics, physics and programming with something really impressive. Creating a working simulation, launching and rendering it is a special kind of something from nothing .

I was inspired by the materials of the lecture SIGGRAPH, just like in the simulation of the fluid . If you want to find more comprehensive information about the above, then take a look at the materials of the course SIGGRAPH 2001 “Introduction to physical modeling” . I give a link to the course in 1997, because Pixar seems to have removed version 2001.

I want to thank Maggie Cai for the wonderful illustration of the couple under the umbrella and for the patience with the painstaking selection of colors, while I can not distinguish blue from gray.

And if you are interested, the illustrations were created in Figma .