Откуда берётся приближение: дискретизация и итерация

Урок разбирает два главных приёма, лежащих в основе почти любого численного алгоритма: дискретизацию и итерацию.

Дискретизация — замена непрерывного объекта (функции, области, времени) конечным набором значений в выбранных точках. Итерация — построение последовательности приближений, сходящейся к ответу.

Дискретизация: непрерывное → конечное

Компьютер не умеет работать с «настоящей» непрерывной функцией $f(x)$ на отрезке — у отрезка бесконечно много точек, а памяти конечно. Поэтому первый шаг почти любого метода: выбрать конечную сетку точек и работать только со значениями функции в них. Хотим взять интеграл — разбиваем отрезок на полоски. Хотим найти производную — берём две близкие точки. Хотим решить дифференциальное уравнение — двигаемся по времени маленькими шагами.

Размер шага сетки $h$ — главный параметр. Чем меньше $h$, тем точнее дискретное приближает непрерывное, но тем больше точек надо обработать. Это фундаментальный компромисс: точность против стоимости.

import math

# Дискретизируем sin(x) на [0, pi] и приблизим площадь суммой прямоугольников
def площадь(n):
    a, b = 0.0, math.pi
    h = (b - a) / n
    s = 0.0
    for i in range(n):
        x = a + (i + 0.5) * h   # середина i-го кусочка
        s += math.sin(x) * h
    return s

for n in [2, 4, 8, 16, 64, 256]:
    приближение = площадь(n)
    print(f"n={n:4}: интеграл = {приближение:.8f}  ошибка = {abs(приближение - 2.0):.2e}")

print("\nточное значение интеграла sin на [0,pi] = 2")

Вывод:

n=   2: интеграл = 2.22144147  ошибка = 2.21e-01
n=   4: интеграл = 2.05234431  ошибка = 5.23e-02
n=   8: интеграл = 2.01290909  ошибка = 1.29e-02
n=  16: интеграл = 2.00321638  ошибка = 3.22e-03
n=  64: интеграл = 2.00020081  ошибка = 2.01e-04
n= 256: интеграл = 2.00001255  ошибка = 1.25e-05

точное значение интеграла sin на [0,pi] = 2

Заметьте закономерность: при удвоении $n$ ошибка падает примерно вчетверо. Это и есть «порядок точности» метода — позже мы научимся его измерять и предсказывать. А пока зафиксируем: дискретизация работает, и измельчая сетку, мы управляем точностью.

Итерация: приближение → лучшее приближение

Второй приём нужен, когда ответ нельзя получить прямой формулой даже после дискретизации. Тогда строят итерационный процесс: берут грубое начальное приближение $x_0$ и по некоторому правилу получают $x_1$, затем $x_2$, и так далее. Если правило выбрано удачно, последовательность сходится к искомому значению.

Классический пример — вычисление квадратного корня методом, который знали ещё в Вавилоне. Чтобы найти $\sqrt{S}$, берём любое $x \gt 0$ и повторяем $x \leftarrow \frac{x + S/x}{2}$ — среднее между $x$ и $S/x$. Если $x$ больше корня, то $S/x$ меньше, и их среднее зажимается к истине.

def корень(S, x0=1.0):
    x = x0
    print(f"старт: x = {x}")
    for i in range(1, 7):
        x = (x + S / x) / 2        # вавилонское правило
        print(f"шаг {i}: x = {x:.15f}")
    return x

import math
S = 612.0
ответ = корень(S)
print(f"\nточное math.sqrt = {math.sqrt(S):.15f}")
print(f"совпало знаков: {ответ == math.sqrt(S)}")

Вывод:

старт: x = 1.0
шаг 1: x = 306.500000000000000
шаг 2: x = 154.248368678629703
шаг 3: x = 79.107997864354715
шаг 4: x = 43.422128682151481
шаг 5: x = 28.758162428779126
шаг 6: x = 25.019538536995718

точное math.sqrt = 24.738633753705962
совпало знаков: False

За шесть шагов из абсурдного старта $x = 1$ мы пришли к 25.02 при истинном 24.74. Ещё пара шагов — и совпадут все 15 знаков. Это поразительно быстро: число верных цифр здесь почти удваивается с каждым шагом (квадратичная сходимость — позже разберём, почему).

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

Любой итерационный метод нуждается в двух вещах: правиле перехода $x_{n+1} = g(x_n)$ и критерии остановки. Правило должно быть «сжимающим» — приближать к ответу. Критерий остановки чаще всего такой: остановиться, когда $|x_{n+1} - x_n|$ стало меньше заданного порога $\varepsilon$, или когда выполнено максимальное число итераций (страховка от зацикливания, если метод вдруг расходится).

x = x0
for k in range(max_iter):
    x_new = g(x)
    if abs(x_new - x) < eps:   # достаточно близко — выходим
        return x_new
    x = x_new
raise RuntimeError("не сошлось за max_iter шагов")

Эта рамка повторяется в десятках методов: Ньютон, секущие, простая итерация, Якоби, Зейдель, степенной метод — все они отличаются только видом функции $g$.

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

Бесконечно измельчать сетку. При очень малом $h$ начинают доминировать ошибки округления (об этом — отдельный урок), и точность не растёт, а падает. Есть оптимальный $h$.
Останавливать итерацию только по числу шагов. Лучше комбинировать: и порог по $|x_{n+1} - x_n|$, и лимит итераций как страховка.
Считать, что итерация всегда сходится. Неудачное правило $g$ или плохой старт могут давать расходимость или зацикливание — сходимость надо обосновывать.

Итоги

Дискретизация заменяет непрерывный объект конечным набором значений; шаг $h$ управляет точностью ценой объёма счёта.
Итерация строит сходящуюся к ответу последовательность приближений по правилу $x_{n+1} = g(x_n)$.
Критерий остановки итерации: малое изменение $|x_{n+1} - x_n|$ плюс лимит шагов как страховка.
Эти два приёма — каркас, на который натянуты почти все методы курса.