Об одном подходе к практической оценке погрешностей численных результатов Текст научной статьи по специальности «Математика»

УДК Б I 9.6(07)

В.П. Житников, Н.М. Шерыхалина, С.С. Поречный

ОБ ОДНОМ ПОДХОДЕ К ПРАКТИЧЕСКОЙ ОЦЕНКЕ ПОГРЕШНОСТЕЙ

ЧИСЛЕННЫХ РЕЗУЛЬТАТОВ

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

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

В габл. 1, 2 представлены результаты разных авторов по задаче о солитоне Стокса

(см. цитируемые в таблицах работы). В скобках после чисел указаны погрешности в единицах последнего приведенного разряда

Таблица I Значения числа Фруда, полученные разными авторами

Число Фруда Автор. Название публикации

1,286(5) Longuet-Higgins M.S., Fenton J.D. Proc. Roy. Soc. London, 1974. A 340

1,2909 Fox M.J.H. Ph.D. thesis. Cambridge univ., 1977

1,290906(15) Hunter J.K., Vanden-Broek J.-M. J.F.M. 1983

1,290889(1) Williams J.M. Phil. Trans. Roy. Soc. London, 1981

1,29089053(8) Evans W.A.B., Ford M.J. Proc. Roy. Soc. London, 1996

1,290890455 Шерыхалина Н.М. ВИНИТИ. №2550-B95-дсп. 1995

1,2908904558 Маклаков Д.В. Euro. Journal of Applied Mathematics. Vol. 13. 2002

1,29089045586 Житников B.I I„ Шерыхалина H.M. Физика волновых процессов. 1998. N" 2-3

1,2908904558634 Житников B.I I., Шерыхалина Н.М. Computational Fluid Dynamics fourn. Vol. 10, N»3. 2(H) 1

1,29089045586335 Шерыхалина H.M., 2005

Таблица 2

Значения параметров солитона, полученные разными авторами

Параметр Значения параметров, полученные

Williams (1981) Evans & Ford (1996) Житников, Шерыхалина Житников, Шерыхалина

вычисленные экстраполированные

Число Ег 1.290889(1) 1,29089053 (8) 1,290890455863 1,2908904558634

Амплитуда 0,833197(2) 0.833199179(94) 0,833199084519 0,8331990845196

Масса 1.970319(2) 1,97032019(47) 1,9703206602 1,9703206601319

Импульс 2,543463 (5) 2,54346767 (46) 2,5434681352 2,5434681351545

Циркуляция 1,714569 1,71456873 (50) 1,7145692406 1,7145692405337

Кинет, энергия 0,535005 (4) 0,535008913(77) 0,5350088360 0,5350088359709

Пот. энергия 0,437670 (3) 0,437672702 (9) 0,43767269344 0,4376726934439

Поли, энергия 0,972675 (6) 0,972681614(85) 0,9726815294 0,9726815294148

(там, где погрешность превышает эту единицу).

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

В связи с этим разработка надежных методов оценки погрешности численных методов весьма актуальна.

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

Z„ = Z + <Ч/",(л) + с/2(я) + ...

... + C[fL(n) + А(п), и;

где z— точное значение: z„ — приближенный результат, полученный при числе узловых точек, равном я;/¡, ...,// — некоторые функции числа узлов.

Для разностных формул численного дифференцирования, квадратурных формул Ньютона — Котеса, разностных методов решения задач для уравнений математической физики и многих других формул и численных методов fj(n) = п ~к>, где kj — произвольные действительные числа. Некоторым другим численным методам, например методу простых итераций. соответствуют функции fj{n) = X~J".

Ест и рассматриваются квадратурные формулы на сетке с постоянным шагом или разностные формулы численного дифференцирования, то в качестве параметра дискретизации вместо я может быть использован шаг сетки h = (b — а)/п (полагая, что нумерация узлов начинается с нуля, в качестве я подразумеваем номер последнего узла сетки). Далее используем я как более общий параметр.

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

Пусть имеется конечная последовательность г),'" = -„,/' = 1,...,/ вычисленных результатов. Тогда можно записать систему линейных алгебраических уравнений

Zn I = ^ + с/|(я,) + с/2(я,) + ... ... + с/£(я,) + Д(я,);

г„2 = г + с/|(«2) + с/2(я2) + ... - + с,Мп,) + Д(я2); .................................................................... (2)

г„/= I+ ^(я,) + с/2(я,) + ... ... + с,/,(п,) + Д(я/):

Решение системы (2) представляет собой задачу идентификации математической модели но результатам численного эксперимента. Если считать Д(я,) неизвестными искомыми параметрами наряду с z, сх, ..., С/, то неизвестных в системе (2) всегда больше, чем уравнений, и она имеет бесконечное множество решений, среди которых точное. Применяя методы регуляризации [2—в] можно получить оценку £ этого точного решения I и оценку погрешности. Однако известные методы регуляризации требуют задания некоторой априорной информации не только о неизвестных Д(я(), но и об искомых z, с^ Результаты такой оценки зависят от этой априорной информации. поэтому могут привести к получению ошибочных оценок.

Известные методы оценки погрешности. Возможен другой известный подход [7]. Рассмотрим линейную комбинацию:

/ / /

Xх'г», =гХл'<

(=1 /=| 1=1

/ / (3)

1=1 1=1

Наложим условия полного подавления компонент погрешности:

¿Л- =1, ¿.г,/;(*,)= 0, ..„¿.V,./-(я,) = 0. (4)

(=1 /=1 /=1

Если / = £ + I и система функций ^(п) линейно независима, то система (4) имеет единственное решение. Отметим, что эта задача в отличие от задачи идентификации имеет точное, независимое от погрешностей Л(я() решение. В результате из (3) получается значение

1 I

2„, = Xх<~п, = г + X ), (5)

которое принимается в качестве приближенной оценки точного значения г. Это значение называют квадратурной формулой Ром-берга [7]. Однако оценку погрешности этого значения найти невозможно в связи с отсутствием априорной информации о А(л().

Пусть 0 = «/л7_1 > 1, = (/ Х, т. е. мы устанавливаем некоторый регулярный способ задания п (чаше всего величина О выбирается равной 2). Для этого случая при /¡{п) = п~к' решение системы (4) может быть получено путем последовательного применения экстраполяционной формулы Ричардсона [7].

Результатом этого на каждом шаге экстраполяции является треугольная матрица экстраполированных значений вида

п 7 п _(1) n r<2) "и

10 zw - - -

20 Z20 -20 - -

40 Z40 (!) -40 _(2) -40 -

80 -80 _(2) "80 -80

(6)

Оценка погрешности проводится по правилу Рунге, т. е.

"20

-О)

. О)

А = z - 5 Л = -v"~ z

-20

40

(I) „(2)

-40

40 '

д(2) _ „(2) _ _(3) 80 — -80 ^80

(7)

Этот процесс построения матрицы экстраполированных по Ричардсону значений (6) по строкам и сравнение (для получения оценки погрешности по правилу Рунге) двух последних элементов последней строки называют методом Ромберга.

Недостатком этого метода является отсутствие обоснованности полученных оценок в связи с неопределенностью А(/?,), поскольку оценки (7) представляют собой частные случаи (5).

Этот процесс можно усовершенствовать.

Численной фильтрацией [1, 8, 9] называется последовательное устранение (подавление) компонент погрешности, т. е. определение отфильтрованных последовательностей г)/', / — 1, ..., Для уравнений (2) фильтрация сводится к линейной комбинации

=а,г),М) причем а, и (3, опреде-

ляются из решения системы двух уравнений

а,+ р,= 1; а/Д«,_,) + (З/Дя,) = 0.

Отсюда получаем формулу фильтрации

_(/> _ .(/-О ■ /ОЧ

* +/М-ШМ-1'

Как правило, рассматриваются случаи, когда I > Ь + 1, тогда уравнений в системе (4) больше, чем неизвестных хг и матрица, аналогичная (6), не является треугольной. Матрица заполняется по столбцам с помощью формулы (8) или других фильтров [8], при этом для оценки погрешности сравнивается каждая пара значений Д^' = г'/' -Для наглядности это можно проиллюстрировать на графике в логарифмическом масштабе (рис. I), где представлены результаты обработки данных, полученных при вычислении второй разностной производной

сГ-у[х) = у(х + ь)-2у(х) + у(.г-ь) + ф)

dx1

/г

При этом шаг сетки И — 1 /п. По оси ординат отложены десятичные логарифмы относительных разностей —Ig5, 5 = по оси абсцисс — десятичные логарифмы п. Разности значений двух столбцов при этом представляются системой точек, которую можно условно соединить между собой некоторой кривой или ломаной. Каждая кривая близка к отрезку прямой, если слагаемые в сумме (1) существенно отличаются по величине. Угловой коэффициент каждого отрезка при /}{п) = п~к> (f/n^^/f/n,) = ф приближенно равен соответствующему показателю kj.

На графике изображено несколько линий, тем самым появляется возможность сравнить их взаимное расположение, заметить несоответствия, неправильность их поведения.

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

12 3 4 5 6 0 12 3 4 5 6

Рис. I. Оценки погрешности при вычислении второй производной: - по правилу Рунге; б — сравнение с точным значением. Прямая у = 19 — 2-1

погрешности Д(л), аналогично задаче, рассмотренной в [1]. При этом применение экстраполяции Ричардсона неправомерно, так как она возможна только при преобладании составляющей, имеющий вид степенной функции с определенным показателем.

В диапазоне, где преобладает нерегулярная погрешность, меняющая знак, в качестве оценки погрешности можно использовать разность пар значений. В эк-страполяционной формуле Ричардсона эта разность делится на достаточно большое число 0*> — 1. Поэтому разница приближенного и экстраполированного результатов оказывается малой и оценка по Рунге дает завышенные по точности результаты. С этим связан сдвиг вверх (при увеличении количества фильтраций) участков линий, соответствующих условию преобладания нерегулярной погрешности. Сравнивая рис. 1, о и 1,6 нетрудно заметить, что оценка по Рунге соответствует реальной погрешности в диапазоне, отделенном хотя бы половиной масштабной единицы от уровня нерегулярной погрешности у = 19 — 2-1 %п.

В [I] предложен локальный критерий для принятия решения о достоверности оценки

погрешности. Отношение ги(/' =|Д,„'*|)/Д<„'||

имеет смысл относительной размытости оценки погрешности. В [I] приведено обо-

снование порогового значения /•)■" = 1/3, т. е. при <1/3 оценка принимается,

а при >1/3 отвергается. Следует отметить. что этот критерий использует правило Рунге — сравнение приближенного результата с более точным, полученным по формуле Ричардсона. На участках с преобладанием нерегулярной погрешности Д(п) надежность этого критерия вызывает сомнения.

Сравнение с эталоном. В отличие от правила Рунге все приближенные значения можно сравнивать с одним числом — эталоном г, которое считается наиболее точным. Правило выбора эталона предложено ниже. Разность Ап=гп-1 представляет собой оценку погрешности приближенного значения I. В этом случае результат расчета с оценкой погрешности представляется в виде интервала г = 1±Д„.

Оценки погрешности зависят от правильности выбора эталона г. Чтобы связать погрешность выбора эталона с величиной размытости, используем условие достоверности оценки:

,(0_-

Преобразуем, обозначив =———,

-Л

2-2 (1) _ =-1 п ^ \ г Гп

*■* п ~ |

<1

г\<

1 -г

(Ю)

Величину 7п назовем относительной размытостью оценки эталона. Величина г\, связана с неопределенностью выбора эталона. Оценим влияние изменения эталона на кривые. Выразим 2 через /*'„ из (9) и подставим в выражение

= -18

= -18

Второе слагаемое отражает влияние выбора эталона. Оценивая отличие линии от прямой, полагаем, что смещение на половину масштабной единицы на графике должно быть достаточно заметным. При этом за границы диапазона неопределенности г'п можно принять —3 < г'„ < 0,75. Тогда из (9)

получим |/;,|<1/3. В случае наличия нерегулярной погрешности, "зашумляющей" верхнюю линию, диапазон неопределенности выбора эталона увеличивается, и это должно быть скомпенсировано соответствующим уменьшением . Но поскольку для оценок используется не три значения, как

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

Для увеличения густоты расположения точек можно использовать несколько последовательностей, начинающихся с разных начальных п. Например, если основная последовательность начинается с п = 5 и далее производится удвоение п (0 = 2), то можно построить еще четыре последовательности, начинающиеся с п = 6, 7, 8, 9. Каждая последовательность фильтруется отдельно, а затем результаты фильтрации совмещаются на одном графике и сравниваются с одним эталоном.

На рис. 2 представлены результаты сравнения с эталоном при оценке погрешностей вычисления второй разностной производной. На рис. 2, а эталон выбран предложенным ниже способом. На рис. 2. б эталон искусственно загрублен в 13-м разряде. Как видно из рис. 2, б, при неправильном выборе эталона возникает характерная картина, связанная с ограничением точности на уровне, соответствующем погрешности выбора эталона.

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

а)

5 1(2/1 6

б)

-1?6

16—

5 1ц// 6

Рис. 2. Сравнение с эталоном при вычислении второй производной при различном выборе эталона: а — по результатам оценки с исключением искомого; б — эталон искусственно загрублен в 13-м разряде

Для уменьшения неопределенности предлагается разделить этапы оценки погрешности и определения эталона. На первом этапе проводится фильтрация по формуле

•'л, ''п, '

устраняющая из последовательности г,',0' неизвестное искомое I- Дальнейшая фильтрация по формуле (8) служит оценкой погрешностей, независимой от выбора эталона (рис. 3). Полученная этим способом оценка позволяет выбрать наилучшие с точки зрения минимума погрешности (или комбинации близких по погрешности значений) соотношения z{lll> и / = /0, которые используются на втором этапе фильтрации для определения эталона

-ч

Рис. 3. Результаты оиенки погрешности с исключением искомого по формуле (9)

Тем самым определена формальная процедура (правило) вычисления эталона, что позволяет избежать неопределенности. Таким способом можно оценить и "увидеть" составляющую погрешности А(п).

Методика решения некорректной задачи оценки погрешности следующая:

фильтрация по формуле (II) и многократная фильтрация на основе формулы (8);

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

определение эталона согласно полученным данным путем многократной фильтрации по формуле (8) без применения формулы (11);

оценка погрешности эталона и размытости оценки.

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

Сравнение результатов вычислений и фильтрации с единым эталоном позволяет в отличие от правила Рунге получать оценки погрешности в условиях преобладания нерегулярной составляющей погрешности Л(л).

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

СПИСОК ЛИТЕРАТУРЫ

1. Шерыхалина Н.М., Поречный С.С. Применение методов многокомпонентного анализа для решения некорректных задач // Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. 6(69). 2008. СПб., 2008. С. 89-96.

2. Морозов В.А. Регулярные методы решения некорректно поставленных задач. М.: Наука, 1987. 240 с.

3. Тихонов А.П., Гончарский A.B., Степанов В.В. и др. Численные методы решения некорректных задач. М.: Наука. 1990. 290 с.

4. Тихонов А.Н., Леонов A.C., Ягола А.Г. Нелинейные некорректные задачи. М.: Наука, 1995. 312 с.

5. Федотов A.M. Некорректные задачи со случайными ошибками в данных. Новосибирск: Наука. 1990. 280 с.

6. Зверев Г.Н., Дембннкий С.И. Оценка эффективности геофизических исследований скважин. М.: Недра, 1982. 224 с.

7. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 2004. 636 с.

8. Шерыхалина Н.М. Методы обработки результатов численного эксперимента для увеличения их точности и надежности // Вестник УГАТУ. Сер. Управление, вычислительная техника и информатика. 2007. Т. 9, № 2 (20). С. 127-137.

9. Жипшков В.П., Шерыхалина Н.М. Обоснование методов филырашш результатов численного эксперимента // Вестник УГАТУ. 2007. Т. 9. № 3. С. 71-79.