Применение фильтрации для обработки результатов численного эксперимента Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ

УДК 519.6

Н. М. ШЕРЫХАЛИНА

ПРИМЕНЕНИЕ ФИЛЬТРАЦИИ ДЛЯ ОБРАБОТКИ РЕЗУЛЬТАТОВ ЧИСЛЕННОГО ЭКСПЕРИМЕНТА

Идеи и методы выделения полезной информации из искаженного помехами сигнала путем его фильтрации, идентификации математических моделей информационного сигнала, принятия решений о достоверности выделенной информации применяются в области анализа результатов численного эксперимента для увеличения их точности и надежности. С помощью этих методов проведено исследование результатов работы двух простейших численных методов. Найдены и объяснены некоторые, на первый взгляд, парадоксальные закономерности. Численный эксперимент ; верификация ; уточнение

ВВЕДЕНИЕ

Необходимость применения информационных технологий в области, где, казалось бы, должны преобладать строгие математические подходы, объясняется следующими причинами. Широкое применение математического моделирования в разных прикладных областях привело к бурному развитию численных методов и программ их реализующих, появлению большой массы результатов вычислений. Однако, проведя анализ состояния дел, приходится признать, что достоверность этих результатов часто вызывает сомнения. Многие методы, применяемые для обоснования достоверности и оценки погрешности, не выдерживают критики, поскольку построены не на строгих научных приемах, а на наивных «бытовых» подходах. Можно привести множество примеров опубликованных численных результатов (например, по одной из задач, решенной многими авторами, сравнительный анализ проведен в [1]), которые содержат ошибки в тех разрядах, которые, как утверждают авторы, должны быть верными. Тем самым эти опубликованные данные несут не информацию, а нечто противоположное.

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

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

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

Для анализа решения многих вычислительных задач можно применять методы, сходные с широко развитыми методами обработки информации технического характера. Представляются весьма перспективными методы анализа численных данных путем экстраполяции этих данных по параметрам дискретизации, в качестве которых могут использоваться шаги по переменным, число интервалов разбиения, число слагаемых в сумме и т. п. Известны методы Ричардсона, Рунге, Ромберга, Эйткена [3,4], которые, как показывает практика, иногда работают, а иногда дают неверные оценки. Возникают естественные вопросы: что кроме численных результатов может быть использовано для подтверждения или отвержения полученных оценок, и могут ли результаты расчета сами себя подтвердить, насколько надежными следует считать такие оценки? Результаты численного решения несут информацию как о решаемой математической задаче, так и о методе решения и используемом оборудовании. Можно

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

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

1. МЕТОД ЧИСЛЕННОЙ ФИЛЬТРАЦИИ

Во многих случаях математическую модель погрешности можно представить в виде [5]:

Z„ — Z = С\П kl + С-2П kl +

+ . . . + CLn~kb + A (n), (1.1)

где г — точное значение; — приближен-

ный результат, полученный при числе узловых точек (или числе слагаемых суммы), равном п; ... ,&£ — произвольные действительные числа ( .

В могут входить не вошедшие в сум-

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

нулю при увеличении .

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

ния можно подойти как к задаче интерполяции зависимости от параметра

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

При увеличении в целое число раз, т. е.

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

числе узлов, равном и соответ-

ственно. Составим линейную комбинацию

4V = ajz,n

/3jZV2 —

= (at

Pi)z-

ат

j'n

/3jn2

a-

(1.2)

и потребуем, чтобы суммарный коэффициент при был равен 1, а при (для определенного ^') равен 0. Тогда ау = — {С}к* — 1) \ =

Отсюда получим формулу фильтрации, которая совпадает с Ричардсо-новской экстраполяционной формулой

ZW = Z

ГС 2

'/ТЫ

Qk

1

(1.3)

Проводя последовательно экстраполяцию по всем парам соседних значений и , получим отфильтрованную зависимость, не содержащую члена с

Cl п

(1) -с)-]П

fc;_l

4+in kj+i

ГД ec|lj

с^п кь + (п)

(1.4)

-к,

Заметим, что отфильтрованная последова-

(1)

тельность Zfij содержит на один член меньше, чем исходная. Если она содержит больше одного члена, то ее также можно отфильтровать, устранив степенную составляющую с n^kl. Операции фильтрации можно повторять последовательно для , если ис-

ходная последовательность содержит достаточное количество членов.

В случае rij ф niQ J_1 задача фильтрации сводится к решению системы L+1 линейных уравнений типа (1.1) численным методом.

Главным ограничением для применения рассмотренного метода является наличие неизвестной составляющей погрешности .

Если зависимость от имеет нерегулярную составляющую , обусловленную по-

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

Ащ < a(l-и +

A(L-l) + A(L-l) QkL+1

AU<-U.

(l.5)

Для метода Ромберга, применяемого к последовательности (1.1) (&* = *), произведение таких множителей ограничено числом, приблизительно равным 8, т. е. метод Ромберга является устойчивым к погрешности исходных данных.

Величину погрешности Д(°) можно оценить экспериментально, этому способствует визуализация результатов фильтрации (рис. 1-3).

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

2. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ

Пусть задана сетка ^ = £о + 3^, з = = 0,1,п, Л = 1/п — шаг сетки. Введем обозначения . Рассмо-

трим разностную производную первого порядка

л

./j+1 - fj h '

(2.l)

= - lg 5 =

Перейдем к численному эксперименту. Рассмотрим в качестве тестового примера вычисление производной функции в

точке .

Результаты удобно представлять визуально в виде графика (рис. l,2), где по оси абсцисс даны значения десятичного логарифма п (х = lgn = — lg Л,), по оси

ординат — десятичного логарифма относительной погрешности формулы (2.l)

[cos(t+l/n)-cos f]w+sin t т е точ-sin I, ’ ' '

ности разностной производной, выраженной в количестве точных десятичных знаков (равной количеству информации в децитах).

На рис. l представлены зависимости при расчетах с двойной точностью (Real*8, Фортран MS DOS). Линия, обозначенная цифрой , соответствует погрешности непосредственно вычисленных значений производной, цифрами , и т. д. — погрешности результатов первой, второй и т. д. фильтрации. Кроме того, на рисунки нанесен график линейной функции у = 19 — х.

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

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

мостью |5(п)| и 5о/п, где до = 10м при М=19. Это качественно совпадает с имеющимися оценками [4].

Рис. 1. Переменные типа Real*8

Рис. 2. Переменные типа Real*4

На рис. 2 даны аналогичные зависимости при расчетах с обычной точностью (Кеа1*4). Видно, что уменьшение длины машинного слова сказывается только при < 6. При этом ограничение можно условно провести с помощью прямой , параллельной

линии 0 и отстоящей от нее на 7 разрядов. Это объясняется ограниченностью точности ( 7) длины мантиссы переменной типа

Кеа1*4. Невозможно повысить точность более, чем на 7 разрядов, так как данные не содержат нужной для этого информации.

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

Ответы на вопросы заключаются в следующем. В Фортране М5 DOS функции вычисляются с повышенной точностью ( )

независимо от объявленной точности аргументов и способа вызова функции. С этой же повышенной точностью производятся арифметические действия с результатами вычислений функций, а также к этой точности преобразуются аргументы функций. Для других трансляторов выводы качественно остаются справедливыми (например, на Паскале точность вычисления функций всегда соответствует типу Extended с М = 18-19). В Фортран-трансляторах, написанных под Windows, нет повышения точности до М = 19, ограничение М=16-17 соответствует длине real*8.

Остается вопрос, ответ на который не следует из сказанного выше. Точность величины определяется объявленной точностью переменных, а следовательно, и точность вычисления отношения (2.1) ограничивается тем же. Почему же точность вычисления производной намного превосходит объявленную точность переменных (см. рис. 2, где объявленная точность переменных — Real*4, М = 7,ав эксперименте точность доходит до 13 значащих цифр)?

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

/ (tj + h + Ад) + Aj+i — / (tj) — Aj ,

h + Ah

= {f (tj) + 2^" ^ ^ + ^/l)2 +

Д,+1- A,)/(h+ Д„) = i/"(0 x

х(,! + Д„)+/'(у^ + %±^.

Тем самым сокращается вместе с

погрешностью, в этом и заключается ответ на вопрос. Отметим, что сумма не мо-

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

определяется не объявленной точностью, а той точностью, с которой вычисляются функции (как было указано выше, для Фортрана MS DOS М = 19).

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

Следует также отметить, что, несмотря на оценку (1.5), величина погрешности округле-

ния практически не зависит от числа проведенных экстраполяций. Это связано с ее хаотичностью и завышением оценок, построенных по принципу «наихудшего случая».

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

На рис. 3 показаны результаты оценки погрешности по правилу Рунге для данных, соответствующих рис. 1. Видно, что оценка по Рунге соответствует реальной погрешности в диапазоне отделенным хотя бы одной значащей цифрой от нерегулярной погрешности. В диапазоне, где преобладает нерегулярная погрешность, правило Рунге дает завышенную точность. Это можно объяснить именно тем, что экстраполяция по Ричардсону практически не меняет величины нерегулярной погрешности. Поэтому разница результатов оказывается малой и оценка по Рунге дает завышенные результаты.

О 1 2 3 4 5 Ign 6

Рис. 3. Оценка по правилу Рунге (в сравнении с рис. 1)

О 1 2 3 4 5 lg п 6

Рис. 4. Оценка по правилу Рунге результатов вычисления второй производной

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

3. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

В качестве примера вычислим методом прямоугольников интеграл

Рис. 5. Численное интегрирование методом прямоугольников

■к/2

sin tdx и h 2_^ / (^i), ti = ih — —.

г=1

(3.1)

7Г

При вычислении суммы накоплением возникает ситуация, описанная в [6] , когда складываются слагаемые разного порядка. Общая оценка погрешности квадратурной формулы на основе модели предельной погрешности [4]

i=l ,2/С2)

— Am + Ar + Aj. <

-M

<co(b - a)h f^L + (b^ a)/max • Ю' f (b-a)2

i J max -lu 5

/max = max \f(t)\, /Wx = max .f{k}(t)

i£[a,6] i£[a,6]

(3.2)

Перейдем к численному эксперименту. Методы фильтрации позволяют устранить погрешность метода численного интегрирования (замены интеграла конечной суммой) с высокой точностью. На рис. 5 показаны результаты численного интегрирования по формуле (3.1). Для сравнения на рисунок нанесена прямая . Тем самым из ре-

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

Теперь исследуем, как накапливается погрешность округления при интегрировании константы (рис. 6, линии 1, 2, 3, 4 соответствуют переменным типа single (real*4), real, double (real*8), extended). Типы переменных соответствуют языку Pascal, в скобках указаны эквивалентные переменные Фортрана. Для сравнения на рисунок нанесен график прямой y = 17 — x. Из сравнения графиков видно, что в этом случае более близкой к реальности оказывается модель предельной погрешности (3.2).

Рис. 6. Численное интегрирование константы. Прямая: у = 17 — х

Рассмотрим интегрирование функции, представляющей собой сумму постоянной и переменной составляющих + A sin ж. На рис. 7,8,9 приведены графики относительной погрешности, соответствующие А = 10-5, Ю-12, Ю-15. Видно, что для величина погрешности округления близка к зависимости , при

уменьшении вначале (для больших погрешность округления приближается к зависимости у=16,5-ж.

0 2 4 6 lgn 8

Рис. 7. Численное интегрирование функции

f(x) = 0,1 + A sin * на отрезке [0, тт/2]. А = 10-5

■Igfi

0 2 4 6 lg п 8

Рис. 8. Численное интегрирование функции на отрезке .

■igs

6 -I-------1--------1---------1--------

0 2 4 0 lg” В

Рис. 9. Численное интегрирование функции на отрезке .

Тем самым из анализа результатов эксперимента можно сделать следующие выводы.

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

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

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

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

Для объяснения этих эффектов рассмотрим в качестве примера процесс суммирования константы, значение которой близко к в десятичной системе счисления, т. е.

С = 0,1с-2 ...См-

В результате суммирования получается величина

к

Sk и с = кс. Sk = 0, 040,2 ... а,м • ltP.

г=1

Отметим, что числа хранятся в памяти ЭВМ в нормализованной форме, т. е. величина порядка выбирается так, что перед запятой должен быть ноль, после запятой — не ноль ( ). Это максимизирует количество

информации, которое может быть сохранено в машинном слове ограниченного размера.

При суммировании может возникнуть погрешность округления, вызванная сдвигом мантиссы меньшего слагаемого при выравнивании порядков.

В табл. 1 даны результаты, характеризующие процесс суммирования: fc — количество слагаемых; Sfc — сумма; Докр — десятичный порядок погрешности округления, возникающей при каждой операции суммирования;

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

Отметим, что при суммировании константы погрешность округления повторяется

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

Т аблица Результаты суммирования константы

к sk V ^окр т

1... 9 ОД... 0,9 0 0 9 0

10...99 1... 9,9 1 Ю-м 90 1Q-M+2

100...999 10... 99,9 2 1Q-M + 1 900 Ю — j\/f _|_ 4

1000...9999 100... 999,9 3 1Q-M+2 9000 Ю -м+6

Тем самым при интегрировании константы погрешность округления определяется в основном суммой последних, больших по порядку одинаковых погрешностей, поэтому и накапливается пропорционально .

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

ВЫВОДЫ

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

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

математике «наихудший случай») фактической реальности. Во вторых, это связано с особенностями программирования на языках высокого уровня (типа Фортран, Паскаль) и использованием ресурсов ЭВМ трансляторами с этих языков.

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

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

1. Sherykhalina, N. M. Application of extrapolation methods of numerical results for improvement of hydrodynamics problem solution / N. M. Sherykhalina, V. P. Zhitnikov // Computational Fluid Dynamics J. 2002. V. 11, No 2. P. 155-160.

2. Кулиш, У. Достоверные вычисления. Базовые численные методы : пер. с англ. / У. Кулиш, Д. Рац, Р. Хаммер, М. Хокс. М.-Ижевск : РХД, 2005. 495 с.

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

4. Волков, Е. А. Численные методы / Е. А. Волков. М.: Наука, 1982. 256 с.

5. Шерыхалина, Н. М. Численная фильтрация данных, искаженных нерегулярной погрешностью / Н. М. Шерыхалина, А. А. Ошмарин // Вестник УГАТУ. 2006. Т.8, № 1 (17). С. 138141.

6. Житников, В. П. Линейные некорректные задачи. Верификация численных результатов / В. П. Житников, Н. М. Шерыхалина, А. Р. Ураков. Уфа : УГАТУ, 2002. 91 с.

ОБ АВТОРЕ

Шерыхалина Наталия Михайловна, доц. каф. компьют. мат. Дипл. инж. (УГАТУ, 1993). Канд. физ.-мат. наук (БГУ, 1996). Иссл. в обл. волновых течений жидкости, уединенных волн, методов оценки погрешности численных результатов.