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

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

УДК 519.6

Н. М. ШЕРЫХАЛИНА, С. С. ПОРЕЧНЫЙ

МОДЕЛИРОВАНИЕ ПОГРЕШНОСТИ И ЧИСЛЕННАЯ ФИЛЬТРАЦИЯ ПРИ РЕШЕНИИ СМЕШАННЫХ ЗАДАЧ ДЛЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

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

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

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

1. СМЕШАННАЯ ЗАДАЧА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ

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

д'2и,

дх2

ди

~dt

начальному условию

и (х, 0) = / (х)

и краевым условиям

и (0, t) = <pi (t) , и (ht) = <Р2 (t) , (t > ç0).

(1)

(2)

(3)

Задача решается методом конечных разностей. На плоскости (ж, t) строится сетка с шагом Л по переменной ж (xi = ih, i = 0,... , п, h = 1/п) и с шагом г по переменной i (tj = jr, ,7=0,... , т, t = Ь/т). Вводятся обозначения .

При решении задачи по явной разностной схеме частные производные в уравнении (1) заменяют разностными аналогами следующим образом

Щ+ij 2?ijj ii-i—i.j îij,j+i

h2

(4)

Обозначим Л = т/К2 (число Куранта). Тогда разностное уравнение запишется в виде

i/.y.j у ] -------- /\ I /.у----- ] ) f 1 2Л) (1>

1-J

Хи.

І+1J ;

і = 1.

,п

1; i = ü, —

(5)

Контактная информация: (347)273-32-00

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

Однако явная схема является устойчивой только при А < 1/2 [3,4]. В противном случае развивается экспоненциальный рост погрешности. Поэтому при использовании явной схемы вычисления приходится вести с шагом по времени, зависящим от .

При решении задачи по неявной схеме в (4) изменяется только вид частной разностной производной по времени

Щ+1^ ^Ч,] Щ—1,] — 1 (6)

Л? _ т ’ (6)

а разностное уравнение принимает вид

Ху, ] (1 1 2А) 11/I \а, \ ] ] •

?' = 1...п — 1: 7 = 0....т .

(7)

В соответствии с (2) и (3) значения = / (ж*), Щ (^), ип^ = 9?2(^'), * =

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

( . Затем, подставляя в (7) , реша-

ем систему уравнений относительно «¿,2 и т. д. для всех . Поскольку матрица си-

стемы уравнений (7) является трехдиагональной, эту систему можно решить методом прогонки.

Применение неявной схемы требует решения системы уравнений, но данная схема является безусловно устойчивой.

2. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

В качестве примера на рис. 1 и 2 приводятся результаты решения задачи для = шта; = /(0), ^2{Ь) = /(1). Время окончания счета состовляет (при

этом отличие решения от предельного при £ —> оо составляет примерно 30%). Погрешности оценивались в точке .

Поскольку искомая функция зави-

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

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

ствующих разным , а затем фильтрации полученной последовательности по .

Результаты фильтрации по приведены на рис. 1 в виде зависимости десятичного логарифма относительной погрешности — <5 от \grri для двух значений п: 10

(рис. 1, а) и 320 (рис. 1, б). Точность, полученная при фильтрации, ограничивается накоплением погрешности округления (прямая у = 16 — | ^т). Эта статистическая зависимость характерна также для методов численного интегрирования функций [1]. При п = (рис. 1, б) ограничение по числу Куранта приводит к ограничению сверху шага по времени (ограничении снизу числа ). Это требует дополнительных затрат ресурсов. В результате максимальная полученная точность, оцененная с размытостью (отношением оценки погрешности оценки погрешности к оценке погрешности [6]) порядка 0,1 ... 0,3, достигает 13 значащих цифр.

о ^-------------т-------------т------------

1 2 3 \grri 4

а

-185

14 12 10 8 6 4 2 0

3,5 4,5 5,5 \%ю 6,5

б

Рис. 1. Зависимость значений погрешностей решения смешанной задачи от шага по времени

\grri = 1 ¿Ь/т: а — п = 10, б — п = 320. Явная схема. Пунктирная прямая у = 16 — 0,5 ^ т

Фильтрация полученной последовательности экстраполированных по т значений в сравнении с эталоном, найденным четвертой фильтрацией (рис. 2, а), показывает завышенные оценки (линия 3, рис. 2, а). Это следует из

рис. 2, б, где приведены результаты сравнения с точным решением

и (х, £) = еГЖ 1 8т7Г.Т .

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

1 1,5 2 1 %п 2,5

а

-1§8

О -I--------------------------------------1----------------------------------1----------------------------------

1 1.5 2 2.5

б

Рис. 2. Зависимость значений погрешностей решения смешанной задачи от шага по х: а — сравнение с эталоном; б — сравнение с точным значением. Явная схема

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

В отличие от рис. 1, диапазон изменения т одинаков для всех , поскольку явная схема безусловно устойчива.

Еще одним отличием является заметное ограничение точности (наличие порога) при возрастании п (рис. 3, в, г). Зависимость этого порогового значения точности от (кривая на рис. 4, б) аппроксимируется прямой = 16,5-2-^тг,т. е. 6 ~ п2. Квадратичнаязави-симость обусловлена погрешностью округления значений при вычислении второй разностной производной в (6), что согласуется с

выводами [1]. Отметим, что для явной схемы этот вид погрешности не имеет места.

-18г

14 12 10 8 6 4 2 0

1 3 5 \gtll 7

-18

и 12 и

1 3 5 \gril 7

б

-18 514 -12 -10 -8 -6 -4 -2 -0 -

1 3 5 І£ПІ 7

в

-І8 5

14 12 10 8 6 4 2 0

1 3 5 \%П1 7

г

Рис. 3. Зависимость значений погрешностей решения смешанной задачи от шага по времени (^ т = ^ Ь/т): а — п = 10, б — п = 80; в — п = = 160, г — п = 640. Неявная схема. Пунктирная прямая

в

а

б

Рис. 4. Зависимость погрешностей решения задачи от шага по х: а — сравнение с эталоном; б — сравнение с точным значением. Неявная схема. Пунктирная прямая у = 16,5 — 2 • ^ п

б

Рис. 5. Зависимость погрешностей решения задачи от шагов по х и і (п = т): а — фильтрация методом Ромберга; б — фильтрация методом Нэвилла. Неявная схема. Пунктирная прямая

а

а

Для уменьшения объема вычислений по получению данных для фильтрации можно вместо поочередного изменения чисел и применить синхронное. На рис. 5 приведены результаты такого исследования. Сравнивая рис. 5, а с рис. 4, а, б, следует отметить, что первый способ позволяет добиться несколько большей точности при меньших п. Уменьшить объем вычислений при синхронном изменении и можно, если применить метод Нэвилла [5], который требует увеличения числа узлов не в два раза, а на константу. Однако этот метод неустойчив к накоплению погрешности округления при проведении фильтраций. Этим и объясняется резкое увеличение погрешности при сравнительно небольших (рис. 5, б).

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

) определяется номером старшей производной, входящей в уравнение (в данном случае ).

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

= в + щ) суммированием пар слагаемых одного порядка. Суммирование организуется по схеме бинарного дерева: складываются пары слагаемых, затем пары пар и т. д. (рис. 6, а).

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

Для применения этого алгоритма при решении смешанной задачи введем новые переменные и преобразуем (5)

+ А (''/ I./ — + '«¿+1,.? + «г-1,0 — 2«г,о + '<4+1,0) •

Результаты решения задачи этим способом приведены на рис. 7.

0

1 2

3

4

а,

а.

а,+ +а

б

Рис. 6. Схема процесса попарного суммирования

-1е5

14

12

10

7

1

" 0

1,5

1 %п

2,5

б

Рис. 7. Результаты применения попарного суммирования для решения смешанной задачи: а — зависимость погрешности от шага по х; б — зависимость погрешности от шага по ^ (сравнение с точным значением). Явная схема

Видно, что накопления погрешности округления не происходит. Погрешность остается на уровне, характерном для погрешности исходных данных 10-15.

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

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

3. АНАЛИЗ РЕЗУЛЬТАТОВ ЧИСЛЕННОГО ЭКСПЕРИМЕНТА

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

Матрица системы уравнений (7), решаемой на каждом временном шаге, имеет меру обусловленности V > 1 [3,4]. Но в этом случае может возникнуть экспоненциальный рост погрешности при повторных шагах по времени 8 ~ е^М1Ут. Этого не наблюдается в реальном вычислительном процессе (см. рис. 3, 4). Тем самым, теория может дать верхние оценки погрешности, которые невозможно использовать для практики.

а

а

Все сказанное выше заставляет искать причину погрешности в другом направлении. В [1] подобная квадратичная зависимость погрешности округления от п наблюдалась при фильтрации результатов вычисления второй разностной производной. При решении смешанной задачи также вычисляется вторая конечно-разностная производная. При этом погрешность округления определяется операцией вычитания близких значений функции "¡j А (в) < ^ (где Д0 - ве-

личина порядка единицы последнего разряда числа в машинном слове), при вычислении второй производной по ж, поэтому увеличивается как п2 при увеличении п. Эта оценка и подтверждается численным экспериментом.

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

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

В явной схеме накопление происходит по другим законам и от порядка производной не зависит.

ВЫВОДЫ

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

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

При решении по явной схеме погрешность определяется механизмом накопления погрешности округления при вычислении сумм [1] (зависит от числа разбиений т по

статистическому закону как у/т) и это накопление устраняется попарным суммированием.

При решении методом конечных разностей численная погрешность определяется погрешностью численного дифференцирования, связанной с округлением значений искомой функции, которая для второй производной оценивается как 10-Мп2 [1], гдеМ — длина мантиссы в десятичных разрядах.

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

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

1. Шерыхалина, Н. М. Применение фильтрации для обработки результатов численного эксперимента / Н. М. Шерыхалина // Вестник УГАТУ (сер. «Управление, вычислительная техника и информатика»). 2007. Т. 9, № 7(25). С. 90-96.

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

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

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

5. Шерыхалина, Н. М. Методы обработки результатов численного эксперимента для увеличения их точности и надежности / Н. М. Шерыхалина // Вестник УГАТУ

(сер. «Управление, вычислительная техника и информатика»). 2007. Т. 9, № 2(20). С. 127-137.

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

7. Болотовский, Ю. И. Опыт моделирования систем силовой электроники в среде OrCAD 9.2 / Ю. И. Болотовский, Г. И. Таназлы // Силовая электроника. 2004. № 1. С. 90-95.

ОБ АВТОРАХ

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

Поречный Сергей Сергеевич, асп. той же каф. Дипл. инж. и магистр по инфор. и прогр. обесп. САПР (УГАТУ, 2006). Готовит дис. в обл. мат. моделирования физ. процессов.