Дискретизация и численные ошибки

Дискретизация и ошибки: почему шаг dt решает судьбу симуляции.

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

Производная через конечную разность

Скорость — это производная координаты: насколько быстро меняется x. На бумаге это предел при бесконечно малом шаге. В компьютере бесконечно малого нет, поэтому берём конечный шаг dt и пишем приближённо: v ≈ (x(t+dt) - x(t)) / dt. Чем меньше dt, тем точнее приближение, но тем больше шагов нужно сделать.

Любую производную можно так заменить. Вторую производную (ускорение, кривизну) приближают через три соседние точки: f'' ≈ (f(x+h) - 2f(x) + f(x-h)) / h². Эта формула — рабочая лошадка всего раздела про волны и поля.

Два вида ошибок

Ошибка усечения (truncation) возникает оттого, что мы оборвали ряд Тейлора: заменили кривую её куском прямой на шаге dt. Она уменьшается, когда уменьшаем шаг. Ошибка округления (round-off) возникает оттого, что компьютер хранит числа с конечной точностью (примерно 15–16 значащих цифр). Она, наоборот, накапливается, когда шагов становится слишком много.

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

import math
def orbit_energy_drift(dt, method):
    GM = 4*math.pi**2
    x, y = 1.0, 0.0
    vx, vy = 0.0, 2*math.pi
    def acc(x,y):
        r=math.hypot(x,y); a=-GM/r**3; return a*x,a*y
    def E():
        r=math.hypot(x,y); return 0.5*(vx*vx+vy*vy)-GM/r
    E0 = E()
    ax, ay = acc(x,y)
    n = int(2.0/dt)
    for _ in range(n):
        if method=="euler":
            ax,ay=acc(x,y)
            x+=vx*dt; y+=vy*dt
            vx+=ax*dt; vy+=ay*dt
        else:
            x+=vx*dt+0.5*ax*dt*dt
            y+=vy*dt+0.5*ay*dt*dt
            ax2,ay2=acc(x,y)
            vx+=0.5*(ax+ax2)*dt; vy+=0.5*(ay+ay2)*dt
            ax,ay=ax2,ay2
    return abs((E()-E0)/E0)
print("Дрейф энергии за 2 года орбиты при разных dt:")
print("  dt       Эйлер         Верле")
for dt in (0.02, 0.01, 0.005):
    e = orbit_energy_drift(dt, "euler")
    v = orbit_energy_drift(dt, "verlet")
    print(f"{dt:.3f}    {e:8.4f}      {v:.2e}")
print("Эйлер копит ошибку; Верле — крошечная и не растёт.")

Вывод:

Дрейф энергии за 2 года орбиты при разных dt:
  dt       Эйлер         Верле
0.020      0.5424      4.51e-09
0.010      0.4425      1.69e-11
0.005      0.3423      6.35e-14
Эйлер копит ошибку; Верле — крошечная и не растёт.

Порядок метода

Метод имеет порядок p, если при уменьшении шага вдвое ошибка падает в 2^p раз. Эйлер — метод первого порядка: половинный шаг даёт примерно половинную ошибку (в выводе 0.54 → 0.44 → 0.34, грубо линейно). Верле — второго порядка по координате, и для орбиты его дрейф ничтожен (10⁻⁹ и меньше). Порядок — это «во сколько раз окупается» уменьшение шага.

Устойчивость

Кроме точности есть устойчивость. Для волн и теплопроводности существует критический шаг: если dt больше его, симуляция не просто неточна — она взрывается, числа улетают в бесконечность за несколько шагов. Это условие Куранта (CFL), мы встретим его в разделе про волны. Устойчивость и точность — разные требования: метод может быть устойчив, но неточен, и наоборот.

Как работает под капотом

Под капотом каждый шаг — это разложение Тейлора, оборванное на каком-то члене. Эйлер удерживает только первый член (линейное приближение), Верле — до второго (учитывает кривизну). Чем больше членов удержано, тем выше порядок и точнее шаг, но тем больше вычислений. Поэтому методы высокого порядка (Рунге-Кутта 4-го) делают несколько оценок ускорения за один шаг.

Частые ошибки

Думать, что мелкий шаг всегда лучше. За некоторым пределом растёт округление и время счёта, а точность не улучшается.
Игнорировать устойчивость. Для волн/тепла превышение CFL рушит симуляцию мгновенно, и это не «баг в коде».
Сравнивать методы по одному шагу. Разница видна на длинной дистанции: Эйлер копит ошибку, симплектические методы — нет.

Итоги

Производные заменяют конечными разностями; цена — ошибка дискретизации.
Есть ошибка усечения (падает с dt) и округления (растёт с числом шагов).
Порядок метода показывает, как быстро падает ошибка при уменьшении шага.
Устойчивость — отдельное требование: превышение критического dt рушит симуляцию.