Источники погрешностей: модель, метод и машинное округление

Калькулятор уверенно показывает 0.3, а язык Python — 0.30000000000000004. Кому верить? И второй сюрприз: программа посчитала координату моста с пятнадцатью знаками после запятой, а через месяц пролёт «гулял» на сантиметры. Численные методы начинаются с честного признания: любой расчёт врёт. Вопрос лишь в том, насколько, откуда и можно ли этим враньём управлять.

Три источника погрешности: модель, метод, округление#

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

• Модель: маятник описывают формулой $T=2\pi \sqrt{l/g}$ без трения и без растяжения нити. Формула точна, а реальный маятник от неё отличается — это вина не математики, а упрощений
• Метод: производную заменяем разностью $\frac{f(x+h)-f(x)}{h}$ — даже при идеальной арифметике ответ отличается от $f^{\prime }(x)$ на величину порядка $h$
• Округление: $\frac{1}{3}$ в десятичной записи бесконечно, а в памяти компьютера под него отведено конечное число разрядов — хвост обрезан

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

Больше всего подлости — в первой. Формула пружины по закону Гука $F=-kx$ идеальна, пока деформации малы; вытяните её вдвое дальше предела пропорциональности — и расчётная жёсткость уйдёт в никуда, хотя арифметика будет безупречной. Поэтому опытный расчётчик сначала спрашивает: насколько честна сама формула? И только потом копает в машинные биты. Диагностика, какой из трёх этажей портит ответ, — навык, и собирается он из простых проверок.

Абсолютная и относительная погрешность: значащие цифры#

Пусть $a$ — точное значение, $\tilde{a}$ — приближённое. Абсолютная погрешность — это величина $\Delta =|a-\tilde{a}|$. Относительная погрешность — отношение $\delta =\frac{\Delta }{|a|}$, обычно в процентах; почему именно она судит о качестве приближения независимо от масштаба — словарная статья. Абсолютная погрешность отвечает на вопрос «на сколько врём», относительная — «на сколько процентов».

Пример из лабораторной. Масса детали: $12.0$ кг с точностью $0.05$ кг. Тогда $\Delta =0.05$ кг, а $\delta =\frac{0.05}{12}\approx 0.4\%$. Толщина плёнки: $0.02$ мм с точностью $0.005$ мм — абсолютная погрешность в десять раз меньше, но относительная $\frac{0.005}{0.02}=25\%$: измерение никуда не годится. Абсолютные цифры льстят, относительная — судит.

Отсюда и правило значащих цифр: число $\pi \approx 3.14159$ с пятью значащими цифрами гарантирует относительную погрешность порядка $10^{-5}$, и никак не лучше. Записывая $0.00120$, мы обещаем три значащих цифры (1, 2, 0 — завершающий ноль значим!), то есть ошибку около $10^{-3}$ относительных процентов. Ответ с восемью знаками после запятой, полученный из данных с двумя, — враньё в лучших традициях отчётности.

Машинная арифметика: почему 0.1 + 0.2 ≠ 0.3#

Компьютер хранит дроби по стандарту IEEE 754. Тип double занимает 64 бита: 1 бит на знак, 11 бит на порядок, 52 бита на мантиссу (плюс неявная ведущая единица — итого 53). Это значит: любое число живёт в памяти в виде $\tilde{a}=\pm m\cdot 2^k$, где мантисса $m$ имеет ровно 53 двоичных разряда. Минимальный зазор между соседними числами около единицы — машинное эпсилон $\epsilon =2^{-52}\approx 2.22\cdot 10^{-16}$.

Десятичная дробь $0.1$ в двоичной системе — бесконечная периодическая дробь $0.0\overline{0011}$, как $\frac{1}{3}$ в десятичной. Компьютер обрезает её на 53-м бите. И $0.1$, и $0.2$ хранятся с крошечной ошибкой; при сложении ошибки складываются, и результат оказывается ближайшим к истине машинным числом — но не тем, которое вы ждали:

• $0.1+0.2$ даёт $0.30000000000000004$ — проверьте в консоли браузера
• Десять раз подряд сложить $0.1$: получится $0.9999999999999999$, а не единица
• $0.145\cdot 100=14.499999999999998$ — реальный источник копеечных расхождений в бухгалтерских системах
• $1-0.9=0.09999999999999998$

53 бита мантиссы — это примерно 15–16 верных десятичных цифр. Для чертежа моста за глаза. Но ошибки перемножаются при миллиардах операций, и потому правило одно: числа с плавающей точкой не сравнивают через равенство. Пишут $|a-b|<\epsilon$ с разумным $\epsilon$ — например, $10^{-9}$.

Как погрешности гуляют по арифметике#

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

Числа. Пусть $a=10\pm 0.1$ (то есть $\delta a=1\%$) и $b=20\pm 0.1$ ($\delta b=0.5\%$). Произведение $a\cdot b=200$ с относительной погрешностью $1.5\%$, то есть $200\pm 3$. Частное $b/a=2$ унаследует те же $1.5\%$: $2\pm 0.03$. Ничего страшного.

А теперь вычитание близких чисел: $a=100\pm 0.1$, $b=99.9\pm 0.1$. Разность $a-b=0.1$, но абсолютная погрешность складывается: $0.1+0.1=0.2$. Относительная погрешность разности — $\frac{0.2}{0.1}=200\%$. Два точных числа с точностью до одной тысячной процента дали ответ, который не заслуживает доверия вообще. Умножение и деление бережут относительную точность, вычитание близких чисел — убивает её. Запомните это как правило арбитра: судит не количество знаков, а относительная погрешность результата.

Катастрофическое сокращение: чем опасно √(x+1) − √x#

Классика жанра — функция $f(x)=\sqrt{x+1}-\sqrt{x}$. Математически она безобидна и при больших $x$ ведёт себя как $\frac{1}{2\sqrt{x}}$. Но посчитайте её «в лоб» на double при больших $x$ — и посмотрите, что выдаёт машина:

Что происходит? При $x=10^{16}$ соседние машинные числа отстоят друг от друга на 2. Значит, $x+1$ округляется обратно в $x$, корни совпадают, и разность равна нулю — при том что истинный ответ $5\cdot 10^{-9}$. А при $x=10^{12}$ вычитание $\sqrt{x+1}-\sqrt{x}\approx 5\cdot 10^{-7}$ уничтожает двенадцать старших цифр, и остаются младшие — те самые, что были испорчены округлением корней. Это и есть катастрофическое сокращение (catastrophic cancellation): близкие числа вычитаются, значащие цифры гибнут.

1. 1
   Умножаем и делим на сопряжённое выражение: $\sqrt{x+1}-\sqrt{x}=\frac{(\sqrt{x+1}-\sqrt{x})(\sqrt{x+1}+\sqrt{x})}{\sqrt{x+1}+\sqrt{x}}$
2. 2
   Числитель сворачивается по разности квадратов: $(x+1)-x=1$
3. 3
   Получаем устойчивую формулу $f(x)=\frac{1}{\sqrt{x+1}+\sqrt{x}}$ — без единого вычитания близких чисел
4. 4
   Проверка на $x=10^{16}$: машина возвращает $5.0000000\cdot 10^{-9}$ — все обещанные цифры на месте

Тот же приём спасает квадратное уравнение: когда $b^2\gg 4ac$, корень $x=\frac{-b+\sqrt{b^2-4ac}}{2a}$ считают через $x=\frac{2c}{-b-\sqrt{b^2-4ac}}$, а второй корень — по Виету. Мораль: прежде чем винить компьютер, поищите алгебраическое преобразование, убирающее вычитание близких величин.

Неустойчивые алгоритмы: ошибка, которая растёт#

Даже без вычитания близких чисел алгоритм может сам раздувать ошибку. Возьмём последовательность $y_0=1$, $y_1=\frac{1}{3}$ и рекуррентную формулу $y_{n+1}=\frac{13}{3}{y}_n-\frac{4}{3}{y}_{n-1}$. Точное решение — $y_n=3^{-n}$: подставьте и проверьте, коэффициенты подобраны так, что корни характеристического уравнения равны $\frac{1}{3}$ и $4$. Но начальные данные хранятся в double, и $\frac{1}{3}$ обрезано на 16-й цифре. Эта крошечная примесь включается в общее решение с множителем $4^n$ — и растёт вчетверо каждый шаг:

Формально мы делали всё правильно: четыре арифметических действия, никаких вычитаний близких чисел. Но алгоритм неустойчив: он усиливает входную ошибку, и она растёт экспоненциально — как множитель $4^n$ никуда не денешь. Устойчивая альтернатива здесь — считать $y_n=3^{-n}$ напрямую или идти рекуррентность в обратную сторону, где ошибка делится на 4, а не умножается.

Отсюда практический вывод, ради которого всё затевалось: точность расчёта ограничивает слабейшее звено цепи «модель → метод → арифметика → алгоритм». Красивые 15 знаков от железа не спасут, если метод даёт ошибку $10^{-3}$, а алгоритм усиливает её на каждом шаге. Поэтому в численных методах сначала оценивают погрешность каждого этажа — и только потом верят ответу. Как диагностировать вредителя:

• Уменьшили шаг, а ошибка упала — главная погрешность метода: алгоритм в порядке, вопрос в дискретизации
• Шаг уже микроскопический, а ошибка стоит на месте или растёт — вы упёрлись в округление: меняйте формулу, а не шаг
• Ответ не сходится с экспериментом при идеальной сходимости — переоцените модель: где-то в ней нет существенного фактора
• Ошибка взрывается со временем или номером шага — ищите неустойчивость: где-то в решении спрятан растущий множитель

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