Выбор шага численного интегрирования

Рассмотрим, как сейчас выполняется расчет нового положения системы.

void Rope::update()
{
    for (unsigned i = 0; i < numberOfParticles; i++)
        particle[i].clearAccumulator();

    // Находим длительность последнего кадра, в секундах
    double duration = (double) TimingData::get().lastFrameDuration * 0.001;
    if (duration <= 0.0) return;

    // Выполняем шаг моделирования
    registry.updateForces(duration);

    for (unsigned i = 0; i < numberOfParticles; i++)
        particle[i].integrate(duration);

    Application::update();
}

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

double currentTime = time_in_seconds();

while (!quit) {
    double newTime = time_in_seconds();
    double frameTime = newTime - currentTime;
    currentTime = newTime;

    integrate(state, frameTime);

    render(state);
}

Здесь функция integrate() выполняет численное интегрирование и определяет новое положение системы state, а render() — отображает это состояние на экране. Текущее время измеряется высокоточным счетчиком time_in_seconds() (вроде того, что мы используем в Timing). frameTime — длительность кадра.

Привязка шага интегрирования к длительности кадра приводит к тому, что результаты моделирования будут отличаться для машин с разными видеосистемами. Например, частота обновления экрана 75 Гц (и соответствующий ей шаг интегрирования) может оказаться достаточной для качественных результатов, тогда как при 60 Гц мы получим "разваливающуюся" картинку, вроде той, что получили выше. Подробнее, эти проблемы рассмотрены здесь.

Решение проблемы: фиксированные шаги по времени

Примем, что видеосистема компьютера "производит" время некоторыми порциями (кадрами), длительность которых от нас не зависит. Необходимо, чтобы система моделирования "потребляла" это время фиксированными порциями (шагами интегрирования) заданной длительности, то есть выполняла следующее:

const double dt = 0.001;

double currentTime = time_in_seconds();
double accumulator = 0.0;

while (!quit) {
    double newTime = time_in_seconds();
    double frameTime = newTime - currentTime;
    currentTime = newTime;

    accumulator += frameTime;

    while (accumulator >= dt) {
        integrate(state, dt);
        accumulator -= dt;
    }

    render(state);
}

Обратите внимание, что при таком подходе к выбору шага у нас остается некоторое "недомоделированное" время в конце кадра, меньшее dt. Но это время не пропадает, а переносится на следующий кадр при помощи переменной-аккумулятора (accumulator).

При выборе шага интегрирования мы можем столкнуться с проблемой так называемой "спирали смерти". Она возникает, когда в погоне за точностью задан столь малый шаг интегрирования, что увеличивается время расчета, а, следовательно, длительность кадра, в течение которого расчет выполняется. Увеличение длительности последнего кадра приводит к тому, что в следующий раз для расчета потребуется еще больше времени, и т. д. Для преодоления этой проблемы после определения frameTime поставим проверку:

    if (frameTime > 0.01)
        frameTime = 0.01;

Теперь длительность кадра в любом случае не будет превосходить некоторой заданной величины (0,01 с).

Реализуем рассмотренный подход.

void Rope::update()
{
    // Находим длительность последнего кадра, в секундах
    double duration = (double) TimingData::get().lastFrameDuration * 0.001;
    if (duration <= 0.0) return;

    if (duration > 0.01 )
        duration = 0.01; // максимальная длительность кадра, во избежание "спирали смерти"

    accumulator += duration;

    while (accumulator >= dt) {
        // Обнуляем действующие силы
        for (unsigned i = 0; i < numberOfParticles; i++)
            particle[i].clearAccumulator();

        // Выполняем шаг моделирования
        registry.updateForces(dt);

        for (unsigned i = 0; i < numberOfParticles; i++)
            particle[i].integrate(dt);

        accumulator -= dt;
    }

    Application::update();
}

Код примера полностью

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

Дальнейшее развитие

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

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

const double dt = 0.001;

double currentTime = time_in_seconds();
double accumulator = 0.0;

State previous;
State current;

while (!quit) {
    double newTime = time_in_seconds();
    double frameTime = newTime - currentTime;
    if (frameTime > 0.01)
        frameTime = 0.01;
    currentTime = newTime;

    accumulator += frameTime;

    while (accumulator >= dt) {
        previousState = currentState;
        integrate(currentState, dt);
        accumulator -= dt;
    }

    const double alpha = accumulator / dt;

    State state = currentState*alpha + previousState*(1.0 - alpha);

    render(state);
}

Теперь на экране будет отображаться не последнее вычисленное состояние системы, а смесь этого состояния с предыдущим, полученная с помощью линейной интерполяции с коэффициентом alpha. Коэффициент смешивания состояний определяется "остатками", накопленными в аккумуляторе: alpha = accumulator / dt. Так, если этот остаток равен 0,9 * dt, то состояние системы в основном определяется текущим расчетным состоянием currentState. Если же, напротив, остаток равен 0,2 * dt, то можно сказать, что на 80% состояние системы определяется предыдущим состоянием previousState. Благодаря интерполяции, программа будет отображать состояние системы, отстающее от текущего момента времени приблизительно на время dt, причем это отставание будет одинаковым для всех кадров.

К решению проблем, связанных с необходимостью использовать малые шаги по времени, можно подойти с другой стороны — используя более совершенные методы численного интегрирования, например метод Рунге-Кутты 4-го порядка. Если погрешность метода Эйлера пропорциональна шагу dt, то погрешность метода Рунге-Кутты пропорциональна dt^4. То есть, чтобы добиться точности 0,0001 можно использовать метод Эйлера с шагом 0,0001 или метод Рунге-Кутты с шагом 0,1 ((0,1)^4 = 0,0001).