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

Вычислительные технологии

Том 18, № 1, 2013

Применение многократной фильтрации при численном решении задач методами теории функций комплексного переменного

В. П. Житников, Н. М. Шерыхалина Уфимский государственный авиационный технический университет, Россия e-mail: zhitnik@ugatu.ac.ru, n_sher@mail.ru

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

Ключевые слова: достоверные вычисления, численная фильтрация, комплексные переменные.

Введение

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

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

1. Математическая модель вычислительного процесса

Математическую модель процесса вычисления многими численными методами некоторой величины г представим в виде зависимости

гп = г + сх/х (п) + С2/2 (п) + ... + с/ (п) + А (п), (1)

где г — точное значение; гп — приближённый результат, полученный при числе узловых точек, равном п; /г,...,/ь — некоторые функции числа узлов. Все входящие в (1) константы и функции могут иметь как действительные, так и комплексные значения.

Для разностных формул численного дифференцирования и квадратурных формул Ньютона — Котеса / (п) = п-к, где к^ — действительные числа, что основано на разложении по формулам Тейлора и Эйлера — Маклорена [3]. В разностных методах решения задач для уравнений математической физики наличие нескольких компонент зависимости погрешности от п указанного вида подтверждается вычислительным экспериментом [4]. Некоторым численным методам, в частности рассмотренным в [1], а также ниже в данной работе, соответствуют функции / (п) = Ап, где | < 1.

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

Пусть имеется конечная последовательность гП0) = , I = 1,...,/, вычисленных результатов. Тогда можно записать следующую систему равенств:

= £ + е/ (п:) + С2/2 (п:) + ... + е/ (п:) + А (щ) ,

Zn2 = £ + С1/1 (п2) + С2/2 (п2) + ... + е/ (п2) + А (п2) ,

¿щ = £ + С1/1 (п/) + С2/2 (п/) + ... + е/ (п/) + А (п/). (2)

Если наряду с £ и с1,... , е^ значение А (п^) считать неизвестными искомыми параметрами, то неизвестных в системе (2) всегда будет больше, чем равенств, и она как система уравнений будет иметь бесконечное множество решений, среди которых есть точное.

Применяя методы регуляризации [5-9], можно получить оценки £ этого точного решения £ и погрешности. Однако известные методы регуляризации требуют задания некоторой дополнительной априорной информации о неизвестных А (п^) и, возможно, об искомых е^-. Оценки, полученные такими методами, зависят от этой дополнительной информации, поэтому их вряд ли можно рассматривать как оценки погрешности.

Во избежание некорректности предлагается разделить задачу на две, включающие соответственно идентификацию математической модели по результатам численного эксперимента и тестирование с помощью известных частных решений или других методов. Первая задача состоит не в определении теоретических значений е^-, а в разложении гп на составляющие (компоненты) по известному заранее или определяемому экспериментально базису, т.е. в представлении зависимости гп в виде (1). Хотя математическая запись такой задачи по форме совпадает с (2), однако при этом А (п) и все другие входящие в (1) компоненты приобретают совершенно другой смысл, поскольку заведомо известно, что в А (п) не входят / (п), ] = 1,... , Ь,и константа. Эту задачу приближённо можно решить с помощью фильтрации — получить эталон £ и оценку погрешности А. При этом найденный интервал неопределённости — А, £ + А может не содержать точного значения например, из-за ошибки в программе.

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

В работе [10] получена теоретическая оценка достоверности (доверительной вероятности) совместного результата решения этих двух задач.

2. Численная фильтрация

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

дится к линейной комбинации zn/ = ®j zn. решения системы двух уравнений

(j-1)

+ вZrji-i), причем aj и в определяются из

aj + ^j = ^ ajfj (ni-i) + fafj (ni) = 0

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

(j)

y(j-i)

+

z(j-i) _ z(j-i)

zni znt-1

Rj-1

Rj

fj (ni-1) fj Ы '

Тогда

(j)

(ajRm + ) Cm ) = Cm )

Rj -1

m

j + 1,

В результате фильтрации имеем новую последовательность, не содержащую компоненты fj(n), которая при условии сохранения вида остальных компонент (точнее, первой из оставшихся) может быть подвергнута повторной фильтрации для подавления следующей компоненты. Для сохранения вида компонент достаточно выполнение условия Rj = const. В случае fj (n) = An (при действительном, или комплексном, Aj) из данного условия следует щ -ni-1 = a = const i, Rj = T-a, причем каждая из показатель-

ных

составляющих представляет собой геометрическую прогрессию. Для fj (n) = n

—kj

из того же условия получаем п^/п^ = Q = сопэ^, Я = Qkj, и (3) совпадает с экс-траполяционной формулой Ричардсона. При этом алгоритм многократной фильтрации представляет собой метод Ромберга [3].

n

n

m

3. Оценка погрешности и правило выбора эталона

Для оценки погрешности часто применяется правило Рунге [3], которое заключается в сравнении приближённого и отфильтрованного значений, т. е. оценкой служит второе слагаемое правой части равенства (3). Следует отметить, что это правило даёт верные оценки при условии преобладания конкретного слагаемого суммы (1) (а именно, исключаемого на данном шаге фильтрации) над остальными. Если же преобладает, например, последнее слагаемое (1), которое может иметь хаотичный характер вследствие наличия погрешности округления, то оценка погрешности, полученная по правилу Рунге, может оказаться намного меньше реальной [2]. В связи с этим в [2] предложено сравнивать вычисленные и отфильтрованные значения с единым числом, названным "эталоном". Ошибка в выборе эталона лишена отмеченного недостатка правила Рунге ("кажущегося уточнения"), вызванного зависимостью оценки от конкретной закономерности изменения погрешности [2]. Однако "экспертный" выбор эталона неудобен вследствие необходимости вмешательства эксперта в процесс оценки в интерактивном режиме.

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

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

1 = 1 — , (4)

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

Отметим, что преобразование (4) изменяет компоненты зависимости (1):

^ 1 = ... + [/? (п»-1) — /? (пг)] е?-1) + ... = ... + [1 — Я-1] е?-1)/ (пг-1) + ..., (5)

однако при | >> 1 это изменение незначительно, чем и объясняется "привязка" результата фильтрации (4) к п = п^-1.

Полученная оценка позволяет выбрать лучшие по минимуму погрешности (или комбинации близких по погрешности значений), соотношения п и ] = . Тем самым приходим к задаче минимизации: при некотором постоянном к = 0,1, 2,... необходимо найти минимум А по г и ] при ограничениях

—А < 4? < А,

—А < ? < Д. (6)

Значения п и ] = полученные при решении задачи минимизации, используются для определения эталона £ = гП?о) путём фильтрации последовательности гП0 = по формуле (3). Таким образом, найдена формальная процедура определения эталона £ = £П?о). Далее в совокупности с оценкой (6) вычисляется некоторое множество [£ — А, £ + А — интервал неопределённости.

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

4. Применение метода фильтрации на примере задачи Хеле-Шоу

В работе [12] решена стационарная задача Хеле-Шоу применительно к электрохимической обработке электродом-инструментом (ЭИ) в виде плоскости с криволинейным (полукруглым радиуса R) выступом A'FCGB'. Электрод-инструмент движется вертикально вниз со скоростью Vet (рис. 1). Межэлектродное пространство заполняется электролитом, включается ток, и происходит растворение обрабатываемой поверхности (анода). При достаточно длительном процессе обрабатываемая поверхность AMDNB приобретает стационарную форму. Асимптотическая величина зазора S известна. Напряжённость на бесконечности слева и справа E0 = U/S, где U — разность потенциалов

Рис. 1. Формы областей, соответствующих межэлектродному пространству: а — на физической плоскости; б — на плоскости комплексного потенциала

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

Образом межэлектродного пространства на плоскости комплексного потенциала является полоса (рис. 1, б). Ввиду симметрии этой области рассмотрим ее правую половину. Функция Ш(£), где £ — параметрическая переменная, областью изменения которой является полукольцо (рис. 2, в), определяется конформным отображением. Производная этой функции

=г

2и ^р

п

те

р

,2т-1

—р т= (р2т-1—1)

р

-(2т-1)/2> т-1

£ т-1 — р(2т-1)/2£

Образом границы ЭИ на плоскости годографа напряженности Е

(рис. 2, а)

является некоторая кривая (при этом на физической плоскости форма границы ЭИ будет представлять полуокружность). Аноду, согласно [12], соответствует разрез по дуге окружности радиуса £0/2.

Рассмотрим функцию (рис. 2, б)

ш

ч 2

г 1п V ^

г

2

—т

Рис. 2. Плоскости годографа напряженности (а), изменения функции (8) (б), параметрического переменного £ (в)

Функцию ш (Z) приближённо будем искать в виде [12]

(К п I -1 Z (Z - v)

Ш (z) = х + г Ь —-2

2 vZ — p2

+ i £ Cmpm (p-™Zm — PmZ

m=1

с действительными коэффициентами Cm. Тогда, учитывая (7) и (8), вычислим диффе-

ренциал

dZ = -

4U

n£0V/C \ vZ — Р:

Z (Z — v) £1 cm(p-

gm = l

г ^m_)

1

- 1

x

x

Vp

Z—p

m=1

p

2m-1

(p2m_ — 1)

(p_

_(2m_ 1)/2 z m_1 — p(2m_1)/2Z _m

Z _

dZ.

10)

Конформное отображение Z(Z) находится с помощью численного интегрирования полученного выше выражения.

Задача решается численно методом коллокаций. Уравнение окружности, задающее форму ЭИ |z (eiCT)| = r2, r = R/S, выполняется в конечном числе точек границы области Z = eiCTm, am = nm/n, m = 0,... , n. Полученная таким образом система нелинейных уравнений решается относительно параметров Cm (m < n), p методом Ньютона с регулированием шага. Численное решение задачи существенно усложняется при p > 0.9, что соответствует r > 30. Поскольку такие решения важны для практических приложений, применение методов уточнения решения имеет большой смысл.

Рассмотрим сначала задачу суммирования ряда в (7) (для n слагаемых).

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

носительных погрешностей — lg 5, 5 =

лП?

(точность данных), ось абсцисс — n.

Результаты непосредственных вычислений (кривая 0) и первой и второй фильтрации этой зависимости (кривые 1 и 2) образуют линии, близкие к прямым, если слагаемые в сумме (1) имеют показательный вид и существенно различаются по величине. Жирными линиями на рис. 3 приведены результаты попарного вычитания (см. (4)), тонкими — данные сравнения с выбранным значением эталона. Отсутствие существенных различий в графиках двух этапов, несмотря на некоторое изменение компонент погрешности

Рис. 3. Оценка погрешности суммирования ряда в (7): а — p = 0.9, б — p = 0.95

(5), подтверждает правильность выбора эталона и достоверность оценок. Видно, что с помощью двух шагов фильтрации (исключения двух компонент (1)) точность порядка 14 значащих цифр может быть достигнута при п = 75 ^ 150.

Следует отметить, что для показательных функций /? в (1) при п = п^-1 • 10, г = 2,...,/, формула фильтрации (3) приобретает вид

Л?)

_?-1) _ „(¿-1) г(?-1) I ^ ^¿-1

+ д-1^ _ 1

Для задачи суммирования ряда в (7) в качестве чисел Л? используются произведения р и £ в различных степенях, т.е. Л? в общем случае являются комплексными числами, причём | Л? | < 1.

Для оценки относительной погрешности результатов решения задачи методом кол-локаций рассмотрим значения торцевого зазора зт, вычисленные при п = 5, 6, 7, ... , 60 (рис. 4, обозначения соответствуют рис. 3). В этом случае различие компонент погрешности (5) учтено за счёт видоизменения преобразований (4), (3)

(0)

'(1 — Л1), £

(?)

Л/-1)

+

_0'-1) _ _(?-1) Л--1 —1

1 - Л?

1Л

1, 2,

?+1

В результате коррекции различие кривых, соответствующих двум этапам фильтрации, практически исчезает. Надёжные оценки погрешности эталонов получаются в диапазоне 10-12 ^ 10-13. При этом для г = 32 такая точность достигнута только благодаря фильтрации.

Путём фильтрации выделены две

три компоненты зависимости

причём при

анализе погрешности метода коллокаций основания Л? этих компонент оказались отрицательными.

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

п

3

п

п

Рис. 4. Оценка погрешности параметра зт (г): а — г = 5, б — г = 32

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

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

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

Список литературы

[1] Житников В.П., Шерыхалина Н.М. Уточнение решений сложных вычислительных задач с помощью постпроцессорной обработки численных результатов // Вычисл. технологии. 2008. Т. 13, № 6. С. 61-65.

[2] Житников В.П., Шерыхалина Н.М., Поречный С.С. Об одном подходе к практической оценке погрешностей численных результатов // Научно-техн. ведомости СПбГПУ. 2009. № 3(80). С. 105-110.

[3] Амосов А.А., Дувинский Ю.А., Копченова Н.В. Вычислительные методы. М.: Изд. дом МЭИ, 2008. 672 с.

[4] Житников В.П., Шерыхалина Н.М. Моделирование течений весомой жидкости с применением методов многокомпонентного анализа. Уфа: Гилем, 2009. 336 с.

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

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

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

[8] Федотов А.М. Некорректные задачи со случайными ошибками в данных. Новосибирск: Наука, 1990. 280 с.

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

[10] Житников В.П., Шерыхалина Н.М. Оценка достоверности численных результатов при наличии нескольких методов решения задачи // Вычисл. технологии. 1999. Т. 4, № 6. С. 77-87.

[11] Житников В.П., Шерыхалина Н.М., ПоРЕЧНЫй С.С. Решение задачи идентификации при оценке погрешностей численных результатов // Научно-техн. ведомости СПбГПУ. 2010. № 1(93). С. 105-110.

[12] Житников В.П., МуксимовА Р.Р., Шерыхалина Н.М., Поречный С.С. Численная оценка параметров нестационарных процессов электрохимического формообразования // Тр. ГОСНИТИ. 2010. Т. 106. С. 67-71.

Поступила в 'редакцию 13 марта 2012 г., с доработки — 21 сентября 2012 г.