Точность машинных вычислений. Часть 2.

Введение

Возвращаемся к теме точности машинных вычислений. В первой статье

http://moderncomputingtechnology.blogspot.ru/2014/02/ComputingPrecision.html.

я рассматривал три типа данных, описывающих вещественные числа: float, double, __float128.

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

В данной статье я рассматриваю методы интегрирования систем обыкновенных дифференциальных уравнений (ОДУ), реализованные с помощью типов float, double, long double и __float128. В механике очень часто изучаются колебательные системы, поэтому все тесты я проводил на системах с периодическим решением.

Вид уравнений

Общий вид системы ОДУ первого порядка, разрешенной относительно производных, в векторной форме имеет вид

где x - n-мерная вектор-функция одной переменной, F - заданная n-мерная вектор-функция n+1 переменной. Система, правая часть которой явно зависит от переменной t, называется неавтономной.

Если функция правой части независит от t, то систему называют автономной, и она может быть записана в виде

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

Стоит отметить, что для исследования таких систем разрабатываются различные программные пакеты. Например, в Ижевском институте компьютерных исследований разрабатывается пакет "Компьютерная динамика: хаос" для исследования динамических систем, демо-версию которого можно найти по адресу http://lab.ics.org.ru/lab/page/kompyuternaya-dinamika/.

Методы решения

Для решения ОДУ разработано множество методов. Здесь я рассмотрю четыре явных одношаговых метода: метод Эйлера, метод Эйлера-Коши, метод Рунге-Кутты 4 порядка, метод Вернера.

Метод Эйлера является самым простым и исторически самым первым методом решения ОДУ, был предложен бессмертным механиком, физиком, математиком и астрономом в одном лице Леонардом Эйлером. Метод Эйлера обладает первым порядком точности и записывается следующим образом:

Метод Эйлера-Коши является усовершенствованием предыдущего, обладает вторым порядком точности и состоит из двух этапов:

Метод Рунге-Кутты 4 порядка - наиболее популярный явный метод интегрирования ОДУ и систем ОДУ, записывается в виде

Для многих задач механики достаточно точности метода Рунге-Кутты 4 порядка, но интересно проверить как проявит себя метод более высого порядка.

Метод Вернера (DVERK) относится к семейству методов Рунге-Кутты и обладает 6-м порядком точности. Уравнения метода здесь не приводятся из-за громоздкости. Описание метода и его таблицу коэффициентов можно найти в книге Hairer, Ernst, Syvert Paul Nørsett, and Gerhard Wanner. Solving ordinary differential equations. Vol. 1. Berlin: springer-Verlag, 1991. или попросить коэффициенты у меня.

При описании методов я несколько раз использовал словосочетание "порядок точности". Если точное решение является полиномом степени p, то метод p-го порядка точности теоретически даст точное решение, а решение по методу порядка p-1 будет неточным. К тому же, при уменьшении шага интегрирования, например, в 10 раз, ошибка решения, опять же теоретически, будет уменьшаться в 10^p раз.

Тестирование на автономной системе

Автономная система вида

встречается при исследовании горизонтального движения маятника, если при этом учитывается вращение Земли. Точное решение системы можно найти в книге Э. Камке Справочник по обыкновенным дифференциальным уравнениям, оно имеет вид

Начальные условия были выбраны следующими:

Значениях параметров:

Интегрирование проводилось на отрезке

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

выбор типа данных не влияет результат при шагах интегрирования до 10^-5. Однако, при шаге 10^-6 реализация с типом float начинает терять точность.

Метода Эйлера-Коши (второй порядок точности) ведет себя аналогично методу Эйлера.

Однако float-реализация теряет при шаге 10^-4, использование типов double, long double, __float128 приводит к одинаковым результатам.

Существенные различия в представлении вещественных чисел проявляются при интегрировании системы ОДУ методом Рунге-Кутты 4 порядка

Наименее точный тип float не дает заявленной методом сходимости. Тип double обеспечивает 4 порядок до шага 10^-3 включительно, тип long double - до шага 10^-4 включительно, тип __float128 - для всех представленных шагов интегрирования.

Метода Вернера 6 порядка дает результаты на несколько порядков более точные, чем метод Рунге-Кутты 4 порядка, но и теряет порядок точности быстрее.

float-реализация вообще не дает сеточной сходимости. double- и long double- реализации начинают терять порядок точности с шага 10^-3. __float128-реализация теряет точность только при шаге 10^-6.

Таким образом, видно две закономерности:

1. При фиксированном методе. Чем точнее тип, тем более мелкие шаги интегрирования можно использовать без потери порядка точности.
2. При фиксированном типе данных. Чем выше порядок точности метода, тем больше шаг, при котором начинает теряться точность.

В идеале хотелось бы использовать метод высокого порядка, что требует более точных типов данных, таких как __float128 или хотя бы long double. Но за точность придется платить процессорным временем. На следующем графике представлены затраты процессорного времени на реализацию одного шага интегрирования.

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

Тестирование на неавтономной системе

Теперь рассмотрим неавтономную систему

решение которой задается следующими выражениями:

При вычислении правой части неавтономной системы ОДУ необходимо знать момент времени, для которого производятся вычисления. В процессе счета при постоянном шаге интегрирования момент времени можно определить двумя способами: аналитически, записывая формулу для n-го момента времени,

или рекуррентно

Аналитическое определение момента времени использует операцию умножения, рекуррентное - операцию сложения. С точки зрения математики разницы нет, но в программной реализации выявляются некоторые особенности, которые демонстрируют следующие 9 графиков.

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

Полученный результат позволяют сделать два вывода:

1. При фиксированном методе. Чем точнее тип данных, тем меньше влияет дополнительная погрешность.
2. При фиксированном типе данных. Чем выше порядок точности метода, тем сильнее он чувствует погрешность.

Послесловие

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

О числах с плавающей точкой можно найти краткую информацию по ссылке.

Если Вас заинтересовали реализации методов, пишите: veugene186@gmail.com

Множество кодов по методам решения ОДУ можно найти по ссылке http://www.unige.ch/~hairer/software.html

Автор выражает благодарность к. ф.-м. н. А.О. Казаков за полезные обсуждения и ценные замечания.