Точность машинных вычислений в задачах вычислительной физики

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

Оран Э. , Борис Дж. Численное моделирование реагирующих потоков,

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

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

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

С какой же точность можно представить действительное число в машинной памяти? Для ответа на этот вопрос рассмотрим три типа данных, представляющих вещественные числа.

Тип одинарной точности float занимает 32 бита памяти и имеет относительную точность 7-8 десятичных цифр.
Тип двойной точности double занимает 64 бита памяти и имеет относительную точность 15-16 десятичных цифр. Именно с этим типом программируются методы современной вычислительной механики.
Тип четвертной точности __float128 занимает 128 бит памяти и имеет относительную точность 31-32 десятичных цифр.

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

Вычисление ряда реализовал с помощью библиотеки длинной арифметики GNU Multi-Precision Library или кратко - gmp.

Здесь я опишу процесс установки в ОС Ubuntu Linux 12.04 64-bit.

1. Скачиваем исходники gmp и документацию с официального сайта https://gmplib.org/. После скачивания у меня в папке Downloads появился архив gmp-5.1.3.tar.gz. 

2. Копируем архив во временную папку temp

cp Downloads/gmp-5.1.3.tar.gz temp/

3. Переходим в папку temp и распаковываем архив 

cd ~/temp

tar -xvzf gmp-5.1.3.tar.gz

4. Переходим в папку извлеченную из архива

cd gmp-5.1.3

5. Запускаем скрипт configure

./configure

6 Запускаем утилиту

make

7. Устанавливаем gmp

sudo make install

В папке /usr/local/include появятся заголовочные файлы gmp.

В папке /usr/local/lib появится динамически подключаемая библиотека libgmp.so

8. Обновляем кэш динамических библиотек

sudo ldconfig

Тестовую программу можно скачать по ссылке

https://drive.google.com/file/d/0Byp7G91RaS2aV3c1ZTNRSUVlRWs/edit?usp=sharing

Для ее сборки выполнить в консоли следующее

g++ test1.cpp -O3 -lgmp -lquadmath --openmp -o test1

Тест запускал на стареньком ноутбуке ASUS K50ID. CPU Intel Core2Duo T6670. Запуск программы показал следующие результаты


Float  2.71828 17459 10644 53125
Double 2.71828 18284 59045 09079 55982 98427 64884 23347 47314 45312
Quad   2.71828 18284 59045 23536 02874 71352 66231 43584 21867 19354
GMP    2.71828 18284 59045 23536 02874 71352 66249 77572 47093 69995
Timing:
Float  8.36839e-10
Double 8.76361e-10
Quad   4.30683e-06

Точность вычислений с типом float - 6 знаков после десятичной точки.

Точность вычислений с типом double - 15 знаков после десятичной точки.

Точность вычислений с типом __float128 - 33 знака после десятичной точки.

Затраты процессорного времени на вычисления с типом float и double примерно одинаковы и невысоки. Затраты на вычисления в __float128 на 4 порядка выше. В принципе для инженера может быть достаточно и шести знаков после запятой, для задач вычислительной механики такой точности может уже не хватить.

В начале статьи я говорил о дифференциальных уравнениях, поэтому второй тест я провел на вычислении конечно-разностных аппроксимаций первой производной со вторым порядком точности

и с шестым

В качестве тестовой функции взял показательную

Эталоном считалось значение, вычисленное с четвертной точностью.

Исходник второго теста можно скачать по ссылке

https://drive.google.com/file/d/0Byp7G91RaS2aZ21Ta2VGeENWa1k/edit?usp=sharing

Для его сборки выполнить в консоли следующее

g++ test2.cpp -O3 -lquadmath -o test2

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

Для формулы второго порядка точности получилась следующая зависимость ошибки вычисления производной от шага

Вычисления с типом float с шагом 0.01 дают приемлемую погреность, однако уменьшение шага ведет к росту ошибки.
Вычисления с типом double приводят к росту ошибки при уменьшении шага ниже 1e-05.
В вычислениях с типом __float128 ошибка округления не проявилась.

Для формулы шестого порядка получились следующие результаты

Вычисления с типом float не обеспечивают приемлемой точности.
Вычисления с типом double приводят к росту ошибки при уменьшении шага ниже 0.01.
В вычислениях с типом __float128 ошибка округления начала проявляться при уменьшении шага ниже 1e-05.


Лично для себя я делаю следующие выводы:
1. В численном дифференцировании не использоваться тип float.
2. Использование типа __float128 дает высокую точность, но может потребовать значительных затрат процессорного времени.
3. Буду как и раньше использовать тип double для вычислений, пока в моем распоряжении не появится процессор, выполняющий операции с __float128 за приемлемое время.
4. Для каждой конкретной задачи нужно отслеживать точность вычислений через исследование сеточной сходимости.

В следующей статье о точности вычислений я рассмотрю задачу численного интегрирования обыкновенного дифференциального уравнения.