Симуляция физического мира

http://jamie-wong.com/post/simulating-the-physical-world/
  • Перевод


Как бы вы подошли к симуляции дождя, ну или любого другого продолжительного физического процесса?

Симуляцию, будь это дождь, поток воздуха над крылом самолёта или же падающий по ступенькам слинки (помните игрушку-пружинку радугу из детства?), можно представить, если знать следующее:

  1. Состояние всего в момент начала симуляции.
  2. Как это состояние меняется из одного момента времени в другой?

Под «состоянием всего» я подразумеваю любые переменные данные, которые либо определяют как выглядит окружение, либо же изменения, происходящие с течением времени. Как пример состояния можно привести положение капли дождя, направление ветра, либо же скорость каждой части пружинки слинки.

Если положить, что всё наше состояние окружения это один большой вектор $\vec{y}$, то можно сформулировать нужные нам данные, указанные выше, в следующее:

  1. Найти значение $y_0$, удовлетворяющее $y(t_0) = y_0$?
  2. Найти функцию $f$ такую, что $\frac{dy(t)}{dt} = f(t, y(t))$.

Зачем нам надо хранить состояние всего в одном векторе, я объясню чуть позже. Это один из тех случаев, когда кажется что мы перебарщиваем с обобщением задачи, но я обещаю, в таком подходе есть свои интересности.

Теперь взглянем как можно хранить всю информацию об окружении в одном векторе на простом примере. Допустим, у нас есть 2 объекта в 2D симуляции. Каждый объект определяется своим положением $\vec x$ и скоростью $\vec v$.


Таким образом, чтобы получить вектор $\vec y$, нам надо объединить вместе вектора $\vec x_1, \vec v_1, \vec x_2 ~и~ \vec y_2$ в один, состоящий из 8 компонент, вот так:

$\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}$


Если вас смущает, почему мы хотим найти $f(t, y(t))$, а не начальное определение $\frac{dy(t)}{dt}$, то смысл в том, что производная нам нужна как выражение, зависящее только от текущего состояния, $y(t)$, констант и самого времени. Если же это невозможно, то наверняка какая либо информация о состоянии не была учтена.

Начальное состояние


Чтобы определить начальное состояние симуляции, нужно определить вектор $\vec y_0$. Таким образом, если начальное состояние нашего примера с двумя объектами выглядит приблизительно так:


То в векторном виде это можно представить следующим образом:

$\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}$


Если объединить это всё в один вектор, мы получим нужный нам $\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}$


Производная функция


$\vec y_0$ определяет начальное состояние и теперь всё что нам нужно это найти способ перехода из начального состояния в состояние происходящее мгновение после него, а с полученного состояния ещё чуть дальше во времени и так далее.

Для этого, решим уравнение $\frac{dy(t)}{dt} = f(t, y(t))$ для $f$. Найдём производную от $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}$


Ого! Высокая получилась формула! Но можно привести её в более читаемый вид, если разобьём наш вектор $\vec y$ обратно на составные части.

$\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$ и $\vec x_2$ связаны аналогичными правилами, равно как и $\vec v_1$ с $\vec v_2$, поэтому несмотря на кучу выражений сверху, всё что нам действительно нужно найти, это следующие 2 вещи:

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


От определения этих двух производных и зависит качество симуляции, именно в них вся сила. И чтобы симуляция не походила на программу где всё случается хаотичным образом, можно обратиться к физике за вдохновением.

Кинематика и динамика


Кинематика и динамика — необходимые ингредиенты для создания интересной симуляции. Начнём с самых основ на нашем примере.

За положение в пространстве отвечает $\vec x$, и первая производная положения точки по времени это его скорость $\vec v$. В свою очередь, производная от скорости по времени это ускорение, $\vec a$.

Может показаться, что мы уже нашли нужную нам функцию $f$, т.к. мы уже знаем следующее:

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


И в самом деле мы блестяще справились с $\frac{d \vec x}{dt}$, т.к. $\vec v$ это часть нашего вектора состояния $\vec y(t)$, но нам нужно ещё чуточку разобрать вторую формулу, потому что с $\vec a$ не всё так просто.

Тут нам поможет второй закон Ньютона: $\vec F = m \vec a$. Если предположить что масса наших объектов известна, то переставив переменные в выражении, мы получим:

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


Так, погодите ка, $\vec a$ и $\vec F$ не являются частью $\vec y(t)$, поэтому это сложно назвать продвижением (помните, нам нужна производная функция зависящая только от $\vec y(t)$ и $t$). Но тем не менее мы продвинулись вперёд, потому что мы нашли все полезные формулы которые определяют поведение объектов в нашем физическом мире.

Предположим, что в нашем простом примере, единственной силой, которая воздействует на объекты является гравитационное притяжение. В таком случае, мы можем определить $\vec F$, используя закон всемирного тяготения Ньютона:

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


Где $G$ это гравитационная постоянная $6.67 \times 10^{-11} \frac{Нм^2}{кг^2}$, а $m_1$ и $m_2$ это массы наших объектов (которые, мы предположим, являются константами).

Для создания самой симуляции, также нам понадобится направление и как то указать $r$ через компоненты вектора $\vec y(t)$. Если предположить что $\vec F_1​​$ это сила действующая на первый объект, то можно сделать это следующим образом:

$\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}$


Подытожим. Изменение состояний в нашей системе из двух объектов полностью выражено через переменные $\vec x_1, \vec v_1, \vec x_2, ~и~ \vec v_2$. И изменения такие:

$\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}$


Теперь у нас есть всё, что отличает нашу симуляцию от всех других симуляций: $\vec y_0$ и $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}$


Но как, имея строго определённую симуляцию, превратить её в красивую анимацию?

Если у вас был опыт написания симуляции или игры, то возможно вы предложите нечто такое:

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

Но давайте чуть остановимся и разберём почему это сработает.

Дифференциальные уравнения


Прежде чем мы приступим к реализации, нужно определиться, какая информация у нас уже имеется и что нам нужно. У нас есть значение $y_0$, которое удовлетворяет $y(t_0) = y_0$, так же есть $f$, удовлетворяющее $\frac{dy}{dt}(t) = f(t, y(t))$. А нужна нам функция, которая даст нам состояние системы в любой момент времени. Математически, нам нужна функция $y(t)$.

Имея это ввиду и приглядевшись к $\frac{dy}{dt}(t) = f(t, y(t))$, можно заметить, что это уравнение связывает $y$ со её производной $\frac{dy}{dt}$. Это означает что наше уравнение дифференциальное! Обыкновенное дифференциальное уравнение первого порядка, если быть точнее. Если его решить, то мы найдём функцию $y(t)$.

Задача нахождения $y(t)$ по данным $y_0$ и $f$ называется задачей Коши.

Численное интегрирование


Для некоторых примеров задач Коши можно легко найти ответ аналитическим методом, но в сложных симуляциях аналитический подход может оказаться очень сложным. Поэтому попробуем найти способ поиска аппроксимированного решения задачи.

Для примера возьмём простую задачу Коши.
Дано: $y(0) = 1$ и $\frac{dy}{dt}(t) = f(t, y(t)) = y(t)$. Найти аппроксимированное решение для $y(t)$.

Рассмотрим задачу с геометрической точки зрения и посмотрим на значение и касательную в точке $t = 0$. Из того, что нам дано, имеем $y(0) = 1$ и $\frac{dy}{dt}(t) = y(t) = 1$


Мы пока не знаем как выглядит $y(t)$, но мы знаем что возле точки $t=0$, значение близко к касательной. Теперь постараемся вычислить $y(0 + h)$ для маленького значения $h$, воспользовавшись касательной. Для начала попробуем $h = 0.5$.


Если расписать, то мы приближаем значение $y(h)$ следующим образом:

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


Так, для $h = 0.5, y(h) \approx 1.5 $.
​​
Теперь мы можем продолжить вычислять для других точек. Хотя, конечно, мы нашли не точное значение $y(h)$, но если наше приближённое значение очень близко к точному, то аппроксимированная касательная тоже будет очень близка к действительной!

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



Далее, продвинемся ещё на $h$ единиц вправо по касательной.


Повторим процесс и получим угловой коэффициент касательной $f(t,y(t))=f(1,2.25)=2.25$:


Процедуру можно проводить рекурсивно и для этого выведем формулу:

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


Данный численный метод решения дифференциальных уравнений называется методом Эйлера. Для общего случая шаг x += v * delta_t.

В нашем конкретном случае, пошаговое решение выглядит так:

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


Используя данный метод, результаты удобно представлять в виде таблицы:

$\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}$


Оказывается, у нашей задачи есть красивое аналитическое решение $y(t) = e^{​t} ​​$:

График функции y(t) = e^​t

Как вы думаете, что произойдёт, если в методе Эйлера уменьшить шаг?


Разница между аппроксимированным и точным решениями уменьшается с уменьшением $h$! К тому же, вдобавок к уменьшению шага, можно использовать и другие методы численного интегрирования, которые могут привести к лучшему результату, такие как метод средних прямоугольников, метод Рунге-Кутты и метода Адамса.

Настало время кодить!


С таким же успехом как мы вывели математическое представление описания симуляции, мы можем написать реализацию симуляции программно.

Т.к. я больше всего знаком с JavaScript, и мне нравится ясность, которую добавляют в код аннотации, все примеры будут написаны на TypeScript.

А начнём мы с версии, в которой подразумевали, что $y$ это одномерный массив чисел, прямо как в нашей математической модели.

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 + h f(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)
}

Оперировать с одномерными массивами не всегда удобно, можно абстрагировать функции сложения и умножения процесса симуляции в интерфейс и получить краткую обобщённую реализацию симуляции используя 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 + h f(t_i, y_i)
            const yNext = yi.plus(f(ti, yi).times(h))
            simulationStep(yNext, tNext)
        })
    }
    simulationStep(y0, 0.0)
}

Положительной стороной данного подхода является возможность сконцентрироваться на основе симуляции: что именно эту симуляцию отличает от любой другой. Используем пример симуляции с двумя объектами, упомянутыми выше:

Код симуляция двух объектов
// Состояние симуляции двух объектов в один тик времени
class TwoParticles implements Numeric<TwoParticles> {
    constructor(
        readonly x1: Vec2, readonly v1: Vec2,
        readonly x2: Vec2, readonly v2: Vec2
    ) { }

    plus(other: TwoParticles) {
        return new TwoParticles(
            this.x1.plus(other.x1), this.v1.plus(other.v1),
            this.x2.plus(other.x2), this.v2.plus(other.v2)
        );
    }

    times(scalar: number) {
        return new TwoParticles(
            this.x1.times(scalar), this.v1.times(scalar),
            this.x2.times(scalar), this.v2.times(scalar)
        )
    }
}

// dy/dt (t) = f(t, y(t))
function f(t: number, y: TwoParticles) {
    const { x1, v1, x2, v2 } = y;
    return new TwoParticles(
        // dx1/dt = v1
        v1,
        // dv1/dt = G*m2*(x2-x1)/|x2-x1|^3
        x2.minus(x1).times(G * m2 / Math.pow(x2.minus(x1).length(), 3)),
        // dx2/dt = v2
        v2,
        // dv2/dt = G*m1*(x1-x1)/|x1-x2|^3
        x1.minus(x2).times(G * m1 / Math.pow(x1.minus(x2).length(), 3))
    )
}

// y(0) = y0
const y0 = new TwoParticles(
    /* x1 */ new Vec2(2, 3),
    /* v1 */ new Vec2(1, 0),
    /* x2 */ new Vec2(4, 1),
    /* v2 */ new Vec2(-1, 0)
)

const canvas = document.createElement("canvas")
canvas.width = 400;
canvas.height = 400;
const ctx = canvas.getContext("2d")!;
document.body.appendChild(canvas);

// Текущее состояние симуляции
function render(y: TwoParticles) {
    const { x1, x2 } = y;
    ctx.fillStyle = "white";
    ctx.fillRect(0, 0, 400, 400);

    ctx.fillStyle = "black";
    ctx.beginPath();
    ctx.ellipse(x1.x*50 + 200, x1.y*50 + 200, 15, 15, 0, 0, 2 * Math.PI);
    ctx.fill();

    ctx.fillStyle = "red";
    ctx.beginPath();
    ctx.ellipse(x2.x*50 + 200, x2.y*50 + 200, 30, 30, 0, 0, 2 * Math.PI);
    ctx.fill();
}

// Запускаем!
runSimulation(y0, f, render)


Если подшаманить с числами, то можно получить симуляцию орбиты Луны!

Симуляция орбиты Луны, 1 пикс. = 2500 км. 1 сек. симуляции равна 1 дню на Земле. Пропорция Луны к Земле увеличена в 10 раз

Столкновения и ограничения


Приведённая математическая модель и в самом деле симулирует физический мир, но в некоторых случаях метод численного интегрирования, к сожалению, ломается.

Представьте симуляцию прыгающего на поверхности мячика.

Состояние симуляции можно описать так:

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


Где $x$ это высота мяча над поверхностью, а $v$ его скорость. Если отпустить мяч с высоты 0.8 метра, то получим:

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


Если изобразить график $x(t)$, то получим нечто следующее:


Во время падения мяча производная функции $f$ вычисляется достаточно легко:

$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}$


С ускорением свободного падения, $a = -g = -9.8 \frac{m}{s^2}$.

Но что произойдёт, когда мяч коснётся поверхности? То, что мяч достиг поверхности мы можем узнать по $x=0$. Но при численном интегрировании, в один момент времени мяч может находиться над поверхностью, а уже в следующий под ней: $x_i > 0, x_{i+1} < 0$.

Можно было бы решить эту задачу путём определения момента столкновения $t_{​c} ~ (t_i < t_c < t_{i+1})$. Но даже если этот момент найти, как определить ускорение $ \frac{dv}{dt}$ так, чтобы оно менялось в противоположную сторону.

Можно, конечно, определить столкновение в ограниченном промежутке времени и применить другую силу $F$ на этот отрезок времени $\Delta t$, но гораздо легче определить дискретную константу ограничивающую симуляцию.

А чтобы уменьшить величину проницания мячом поверхности, можно за один тик вычислять сразу несколько шагов симуляции. В совокупности с этим, код нашей симуляции изменится так:


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)
}

И теперь уже можно написать код нашего прыгающего мячика:

Код прыгающего мячика

const g = -9.8; // m / s^2
const r = 0.2; // m

class Ball implements Numeric<Ball> {
    constructor(readonly x: number, readonly v: number) { }
    plus(other: Ball) { return new Ball(this.x + other.x, this.v + other.v) }
    times(scalar: number) { return new Ball(this.x * scalar, this.v * scalar) }
}

function f(t: number, y: Ball) {
    const { x, v } = y
    return new Ball(v, g)
}

function applyConstraints(y: Ball): Ball {
    const { x, v } = y
    if (x <= 0 && v < 0) {
        return new Ball(x, -v)
    }
    return y
}

const y0 = new Ball(
    /* x */ 0.8,
    /* v */ 0
)

function render(y: Ball) {
    ctx.clearRect(0, 0, 400, 400)
    ctx.fillStyle = '#EB5757'
    ctx.beginPath()
    ctx.ellipse(200, 400 - ((y.x + r) * 300), r * 300, r * 300, 0, 0, 2 * Math.PI)
    ctx.fill()
}

runSimulation(y0, f, applyConstraints, 30, render)




Внимание разработчикам!


Хоть у такой модели есть свои плюсы, она не всегда ведёт к производительным симуляциям. По мне, такой фреймворк полезен для представления поведения симуляции, даже если в ней происходит много чего лишнего.

Например, симуляция дождя в начале этой статьи [прим. В оригинальной статье, в начале вставлена красивая интерактивная симуляция дождя, рекомендую посмотреть воочию] не была создана с использованием, описанного в статье, шаблона. Это был эксперимент с использованием Entity–component–system. Исходники симуляции можно найти тут: симуляция дождя на GitHub.

До скорого!


Я нахожу пересечение математики, физики и программирования чем-то действительно впечатляющим. Создание работающей симуляции, её запуск и рендеринг это некий особенный вид чего-то из ничего.

На всё изложенное меня вдохновили материалы лекции SIGGRAPH, точно так же как и в симуляции жидкости. Если хотите найти более исчерпывающую информацию о вышеизложенном, то взгляните на материалы курса SIGGRAPH 2001 «Введение в физическое моделирование». Привожу ссылку на курс 1997 года, т.к. Pixar похоже удалила версию 2001.

Хочу поблагодарить Maggie Cai за чудесную иллюстрацию пары под зонтом и за терпение при кропотливом подборе цветов, в то время как я не могу отличить синее от серого.

А если вас интересует, то иллюстрации были созданы в Figma.
Переводить ли статью «симуляция жидкости» (ссылка дана в заключении статьи)

Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.

Поделиться публикацией
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама
Комментарии 31
  • +1
    Лучше переведите Rigid Body Dynamics motion with constraints, там больше интересного. И метод решения уравнений для статического равновесия. Еще есть интересные статьи как эта математика распараллеливается на gpu.
  • 0

    Реально Круто! Всегда мечтал о таком!


    А реализация через OpenCL или CUDA будет?

    • 0
      Да тут же простенькие симуляции, и четверти ядра хватит :)
    • +3
      У вас луна улетает от земли тихонько.
      • 0
        Я думаю это из-за недостатка Солнца
        • +3
          Нет, всё же это из-за погрешности при численном интегрировании.
          • +4
            Разумеется из-за погрешностей. Метод Эйлера (тем более явный) — решение «в лоб». Для более точного решения стоит использовать хотя бы leap-frog (вычисляем попеременно то координату, то скорость). Будет и энергия сохраняться, и орбита.

            А «большие дяди» используют хотя бы метод Рунге-Кутты 4 порядка. Пишется в 4 строчки, а погрешность — O(h^2).
            • 0
              Вот за этим я всегда прихожу в комменты, за leap-frog первый раз слышу, хотя написал немало методов. Он очевиден и легко «изобретается», но не знал про его свойства.
              Рунге-Кутты 4 порядка всё же имеет 4-й порядок сходимости, соответственно O(h^4). Может вы имели в виду Рунге-Кутты 2-го порядка? Он действительно в 4 строчки.
              • 0
                Да, разумеется, О(h^4). Второпях набирал =)

                Если еще хотите поэкспериментировать с N-body, гуглите в сторону метода потенциалов: вычисляем не силу, действующую на каждое тело, а потенциал (например, гравитационный) в нужной точке. С него потом можно переходить на Practical Mesh и Tree-Code
        • 0
          Тоже хотел про это написать и даже был готов возмущаться, как так можно было написать. А потом подумал, как сделать симуляцию устойчивую к этой проблеме. Не поднять точность, а принципиально избавить метод от подобных ошибок. Конкретно эту задачу решить несложно — задача двух тел решается аналитически. А как быть с тремя телами? Пока могу придумать только количественные решения, например, методы высоких порядков сходимости (Р-К 4).
          • +1
            Консервативные схемы интегрирования нужны. Короче насильно закон сохранения энергии и импульса поддерживать.
            • 0
              Задача трех тел аналитически не решается (строго говоря, нашелся один финн, который смог ее решить при помощи ну уж ооооочень медленно сходящихся рядов, так что это практического смысла не имеет).

              Чуть выше уже написал, что можно делать в случае, когда тел не меньше трех.
          • 0
            Касаемо симуляции прыгающего мячика. Когда-то я моделировал движение турбины, которая могла соприкасаться со статором и изучал хаотические колебания ротора. Модель соударения я моделировал в приближении Герца, то есть в момент удара действует упругая сила. В этом случае, если использовать метод Адамса для решение ОДУ, то решение получается довольно гладким. Вообще для таких задач надо тщательно выбирать метод численного решения.
            • 0

              Интересно, а как все сделано в PowderToy?

              • 0
                Когда писал себе симулятор гравитации обнаружил, что при сближении объектов (столкновение, а не движение по стабильной орбите) все идет в тартары. А именно: объекты после столкновения разлетаются с сильно увеличенной скоростью и нарушением закона сохранения энергии. Ваша реализация точно также ведет себя, если обнулить скорость Луны.
                Есть какие-то способы это обойти?
                • 0

                  Гравитационный манёвр?

                  • 0
                    А как столкновение обрабатывали? Если это было именно столкновение, то уже 100 раз написано в множестве физдвижков и куча статей про это.
                    Если же вы не обрабатывали это отдельно, а это был пролёт очень близко к центру масс, то всё ещё проще. Скорости около центров огромные, а влияние малейших изменений расстояния очень велико, поэтому любые погрешности (как метода, так и вычислительные) резко усиливаются. Выше писали про leap-frog, он не даст супер точности, но с энергией должно быть хорошо.
                    • 0
                      Столкновение никак не обрабатывалось. Строго говоря я не хотел бы вообще их обрабатывать до того момента как пойму, что моя модель корректно работает на гравитации. Про корректность: я ожидал что два идеальных шарика в вакууме (изначально с нулевой скоростью) будут бесконечно мотыляться вдоль одного ограниченного отрезка. За leapfrog посмотрю как оно себя ведет, спасибо.
                      • 0
                        Уточнение: под столкновением в оригинальном комментарии я подразумевал такое сближение объектов, что их координаты при максимально точных рассчетах должны были рано или поздно совпасть. Самого столкновения не было
                        • 0
                          leap-frog поможет на больших расстояниях, Луна перестанет стабильно удаляться. В вашем случае этого может не хватить, потому что рядом с нолём сила уходит в бесконечность и правильно проинтегрировать вряд ли получится. Методы более высоких порядков дадут выше точность, но всё равно ничего не гарантируют. Метод потенциалов, выше, осилит подобное.
                          Но скорее всего в вашей задаче не нужна эта безумная точность в окрестности ноля, для небесной механики достаточно ограничить расстояние радиусом планет. При сближении на такое расстояние посчитать соударение шариков (неупругих с потерями на нагрев).
                          • 0
                            По глупости изначально решил, что метод потенциалов не поможет, но сейчас понял что был не прав. Спасибо.
                            А вот касательно безумной точности вы не совсем правы. Когда клацал по чужим симуляторам (которые «решали» эту проблему через соударение) в некоторых случаях добивался точно такого же нарушения. Либо их писали раки, либо это частичное решение проблемы.
                            Но у меня задача была по моделированию разных вариантов гравитации для собственного любопытства. Тут именно максимальная точность нужна.
                            • 0
                              Либо их писали раки

                              Как бы грубо не звучало, но скорее всего так. В физдвиижках давно научились правильно обрабатывать столкновения. Это же два бильярдных шара, там всё банально считается в лоб. Вывести формулы скоростей после соударения, да тут первокурсник справится! Главное правильно отследить момент удара, а не давать погружаться телам друг в друга за один кадр.
                              Для любопытства это здорово, но планета — не материальная точка, у неё размер есть. Взаимодействия на сравнимых с размером дистанциях уже интереснее. Материальной точке пофиг, а протяжённые объекты получат приливной градиент. Внутри планеты тоже всё будет не так. Ровно в центре гравитация вообще 0, а не бесконечность как в формуле G*M/R^2.
                              Моя мысль в том, что малые расстояния находятся за пределами применимости формул небесной механики, не надо пытаться добиваться точности от них.
                    • 0
                      Видимо, происходит абсолютно упругое столкновение, т. е. не учитывается переход энергии движения в деформацию / разрушение / нагрев. Надо особо это обрабатывать. Даже мячик в симуляции выше не должен бы прыгать вечно.
                      • 0
                        Ответил выше, я некорректно выразился.
                    • +1

                      Автор очень так нехило замахнулся, я погляжу


                      Симуляцию, будь это дождь, поток воздуха над крылом самолёта или же падающий по ступенькам слинки

                      Два из трёх по обозначенной методике моделировать — прямо скажем идея так себе. Дождь вменяемо моделировать только стохастическим процессом, там сплошная теория вероятностей, ибо считать динамику миллиардов капель — занятие для бессмертных. Поток воздуха над крылом самолёта это уравнения Навье-Стокса, которые вообще-то в частных производных и требуют решения смешанной (начально-краевой) задачи, а не задачи Коши. Пружинка — ещё туда-сюда, но по-хорошему там тоже частные производные по пространству вылезают.


                      Схема численного интегрирования тоже выбрана не самая удачная, даже для задачи Коши. Хоть бы в общем виде была разностная схема, а не просто метод Эйлера. Ну и решать большие сильно-обусловленные задачи с высшими производными (которых на практике большинство) в явном виде без контроля погрешности это конечно для больших ценителей, обычно-то берется что-нибудь квазиньютоновское с переменным шагом типа CVODE.


                      В общем какая-то странная статья. Во введении замах на "как численно моделировать всё", а внутри "численное моделирование ньютоновской динамики для самых маленьких". Напоминает студенческий реферат чем-то.

                      • +1
                        Я думаю автор ни на что большее не претендовал, всё же это его статья в его личном блоге, а не публикация в журнале. Для общего представления моделирования, по мне, самое то, как раз для самых маленьких. А введение оно же как рекламный плакат, для привлечения внимания, ничего конкретного он не обещает, опять же, это не научная статья, а его опыт.
                        • +2
                          введение оно же как рекламный плакат

                          Если оно "как рекламный плакат", то это называется false advertising, по-русски "ненадлежащая реклама" (в данном случае недостоверная), нет? Вот я например залез в надежде на универсальный подход, про который может быть чего-то не знаю, а нашёл… ну в общем разочарование. Претензия, безусловно, не к Вам.


                          Для общего представления моделирования

                          Оно, увы, получается совсем не общим, а очень даже узким. Потому что задачей Коши да и даже обыкновенными дифференциальными уравнениями тема не исчерпывается.


                          Вы меня только не поймите неправильно, на уровне "что нового я узнал за лето" статья нормальная. Просто я одно время преподавал численные методы и мне такая подача информации по старой памяти несколько против шерсти, вот я и ною.


                          На правах "выкрика из зала", может, чтобы не получался пересказ Маршака и Пастернака, действительно, взять непосредственно SIGGRAPH'овские конспекты вместо блога "по мотивам"? Подача материала конечно суховата, да и много там букв, но ценность статей кажется больше получится. Вот кстати материалы за 2001-й год, который Pixar никуда вроде не удалил.

                          • 0
                            Если оно «как рекламный плакат», то это называется false advertising, по-русски «ненадлежащая реклама» (в данном случае недостоверная), нет?
                            Возможно, с натяжкой, но в защиту автора скажу, что в начале он ничего конкретного не обещал, только лишь указал на общую модель. У этой статьи, как и у других есть целевая аудитория, не все же помнят или даже проходили численные методы с дифурами. А такая подача как раз и направлена на то, чтобы заинтересовать аудиторию этими темами. И извинюсь за автора, если Ваши ожидания не оправдались :)

                            А за ссылку спасибо. Видимо перенесли материалы курса к себе на страницу. Обязательно передам автору.

                      Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.