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

УДК532.538

ОЦЕНКА ПОГРЕШНОСТИ И ЕЕ ОБОСНОВАНИЕ С ПОМОЩЬЮ ФИЛЬТРАЦИИ ЧИСЛЕННЫХ РЕЗУЛЬТАТОВ, ПОЛУЧЕННЫХ ПРИ РАЗНЫХ ЧИСЛАХ УЗЛОВЫХ

ТОЧЕК СЕТКИ

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

Уфимский государственный авиационный технический университет

Статья поступила в редакцию 15.03.2017

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

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

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

Численная фильтрация. Математическая модель процесса вычисления многими численными методами может быть представлена в виде зависимости

-ki

-k 2

-kr

zn = z + c^n 1 + C2 n 2 +... + cLn L +

A(n)

где точное значение; Zn- приближенный результат, полученный при числе узловых точек, равном п; к - произвольные действительные числа (считается, что к\<к.2<...< кь).

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

Пусть имеется конечная последовательность

i = 1,..., I

z<0> = z

ni n

выгчисленныгх результатов.

Житников Владимир Павлович, доктор физико-математических наук, профессор. E-mail: [email protected] Шерыхалина Наталия Михайловна, доктор технических наук. E-mail: n sher @mail.ru

Соколова Александра Алексеевна, аспирантка. E-mail: [email protected]

Тогда можно записать систему равенств

zn, = z + qn kl + C2ni k2 +... + cLn- kL + A(ni)

"ni

Zn2 = z"

= z + c^^ + C2^k2 + .•• + cLn^kL +A(n2 )

-k

-k

znj = z

= z + cin-kl + C2n- k2 +... + CLn- kL + A(n J )

(2)

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

ныгх последовательностей z(nJ'), ;=1,2,... ( - порядко-

пг

вый номер фильтрации). Для равенств (2) фильтрация сводится к линейной комбинации

(1)

) = a Jj-1) +ß Jj-1)

n I n I J n

'i-1

, причем a j и ß j опреде

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

-k

-k

а з +Р з = 1 а }пгЧ ]п1 ' = 0

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

zС-1) - х^-1) 2(Г) = z О'-1) + пг пг-1

П' П' -!

ni

Q = —— = const ni-1

c(j) = c(j-i) Qkj - Qkm

Qkj -1

, m=j+1,...„L (3)

В результате фильтрации получается новая последовательность, не содержащая компоненты

_

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

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

А

л j)

- z„

. Отношение оценок последовательных

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

Постановка задачи.

Пусть струя идеальной жидкости движется вдоль горизонтальной прямолинейной стенки ADC (рис. 1а). Сила тяжести действует вертикально вниз. Рассматривается решение типа уединенной волны наибольшей амплитуды, при этом на вершине волны образуется излом с внутренним углом 2п/3 (волна Стокса). В точке B скорость жидкости равна нулю. Скорость на бесконечности - V», асимптотическая толщина струи - h. Давление Р на свободной границе равно атмосферному Ро.

а б

Рис. 1. Формы области, соответствующей течению, на плоскостях: а - физической Z = X + iY ; б - комплексного потенциала

На свободной поверхности ABC' значение модуля вектора скорости течения V связано с высотой точки Y уравнением Бернулли при P=Po

2Y

+ —— = const Fr 2 h

Fr =

, ^. (4) где § - ускорение свободного падения; Рг - число Фруда.

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

№(х) = 9 + = 0™{х) = _в(х _ *), (5)

где 9, у - потенциал и функция тока течения; Q = - расход жидкости в струе.

D

С'

А'

а б

Рис. 2. Формы областей, соответствующих течению на плоскостях: а - параметрического переменного х; б - логарифмического годографа скорости

Эта задача решалась многими исследователями [1, 2, 5-13] в течение ряда лет. В табл. 1 представлены численные результаты, полученные разными авторами, по числу Фруда для солитона Стокса. Число Фруда позволяет сравнивать условия волнообразования для разных глубин. В скобках после чисел указаны погрешности в единицах последнего приведенного разряда там, где погрешность превышает эту единицу. Следует отметить, что в цитируемых работах авторы, как правило, декларируют оценку погрешности на уровне единицы последнего приведенного разряда. Сравнение разных результатов показывает, что авторы часто ошибаются в оценке погрешности, причем иногда ошибаются в 10 и более раз. Применение фильтрации позволяет избежать таких ошибок.

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

2

Fr=1.286 (5) Longuet-Higgins M.S., Fenton J.D. Proc. R. S. London. 1974. [5]

Fr=1.2909 Fox M.J.H. Ph.D. thesis. Cambridge univ. 1977. [6]

Fr=1.290906 (15) Hunter J.K., Vanden-Broek J.-M. J.F.M. 1983. [7]

Fr=1.290889 (1) Williams J.M. Phil. Trans. Roy. Soc. London. 1981. [8]

Fr=1.29089053 (7) Evans W.A.B., Ford M.J. Proc. R. Soc. London 1996. [9]

Fr=1.290890455 Шерыхалина Н.М. ВИНИТИ. №2550-В95-деп. 1995 [10]

Fr=1.2908904558 Maklakov D.V. Euro. Journ. of Applied Mathematics. V. 13, 2002. [11]

Fr=1.29089045586 Житников В.П., Шерыхалина Н.М. Физика волновых процессов, № 23, 1998. [12]

Fr=1.2908904558634 Sherykhalina N.M., Zhitnikov V.P. Computational Fluid Dynamics Journ. V. 10, N 3, 2001. [13]

Fr=1.29089045586335 Житников В.П., Шерыхалина Н.М. Уфа: Гилем. 2009. [1]

Fr=1.290890455863341 Zhitnikov V.P., Sherykhalina N.M., Sokolova A.A. Mediterranean Journ. of Soc. Sci. 2015. [2]

Цель работы: подтверждение последнего полученного числа за счет повышения точности решения.

Метод решения задачи. Пусть Z=X+iY, z=Z/h. Решение удобнее получать, если использовать плоскость логарифмического годографа скорости (функции Жуковского, рис. 2б)

dw

га = 1 ln — = 0 + 1т dz

(6)

®(x) = ®i(x) + ®2 (х)

(7)

ГО1 (х) = 01 (х) + Т1 (х) = -i J Re ®1 (ст)

sh пст

ch пст - ch пх

-da

При х ^ + г0, согласно формуле Сохоцко-го [4], получается главное значение (^р.) интеграла

™1(стт ) = -i' v.p. J Re <»1 (ст)

sh пст

ch пст - ch пст.

-da + Re <»1 (стт )

Численное интегрирование производится по двухточечной квадратурной формуле Гаусса 4-го порядка точности. Представим краевое условие (4) в дифференциальной форме

3т dx 1

eJX---- sin 0 = 0

da Fr2

(8)

Для выделения особенности в точке х = 0 определим функцию

f (х) =

1 - е"

(1 + 1е-*х/2 )2

где 9 - угол наклона вектора скорости жидкости к оси абсцисс т=1п(У/У°°).

Функция ю (х) должна удовлетворять следующим условиям:

1) при Reх=0, 0 < 1т х < 1 действительная часть Re ю(х) = 9 = 0;

2) при 1т х = 1 действительная часть Яе ю(х) = 0 ;

3) при 1тх=0 действительная и мнимая части ю (х) связаны уравнением (1);

4) при х = ст + г0, ст ^ ±0 действительная часть Яе ю(х)^ ±п/6;

5) при х ^ ж величина ю(х) ^ 0 .

Задача решается видоизмененным методом Леви-Чивиты [14]. Представим ю(х) в виде суммы

которая имеет действительные положительные значения при х = ст + г и при х = 0 + ги (0<и< 1), стремится к 1 при Reх^■+ж>, стремится к 0 при х^0. Функцию Ю2 (х) представим в виде:

га- (х) = i ln f (х) + iC1 (f (х))Р-1

(9)

Тогда

d ¡ \ п 1 п 1

—ю2 (х) =----+1---+

dX 6ch пх 3 sh пх

2

(

+ipq

1 - е"пх

■чР-1

(1 + 1е _пх/2 )

е~п'х!2 +1 -пх/2 е +1

пх/

(1 + 1е~пх/2 )3

При х = ст +10

га2 ^

(ст +10) = - arctg е 2 + - ln

-пСТ2 + i ь1^

3 1 + е-пст

где Ю2(х) - функция, удовлетворяющая условиям 1, 2, 4, 5.

В результате такого представления га1 (х) — 0 + 1т1 (0) при х — 0, т.е. функция га1(х) является непрерывной на границе. Функция Ю1(х) получается следующим образом. Будем искать на границе х=ст в узловых точках стт (m=0,...,n) значения Re га1(стт )=0m. В силу предыдущего замечания 00 (0) = 0. Кроме того, положим 0n (стп) = 0, поскольку функция 0(ст) экспоненциально убывает при ст — да. В промежуточных между узлами точках функция Re га (ст) получается с помощью интерполяции узловых значений 0т (ст m) кубическим сплайном, дважды непрерывно дифференцируемым. Для восстановления функции ®1 (х) используем формулу Шварца [1]

+ C1

(

Л -пстЛР 1 - е

1 + е"

i(zp arctg е"пст/2 )+

+ 1С

1 - е"

ЛР

1 + е"

cos(2p arctg е пст^2 )- С =

= 02 (ст)+ 1т2 (ст). При х = ст +10, ст —^ 0 + 0 для 0.5 <Р< 1 имеют место оценки: f (ст +10) = ^+ O (ст2),

02 (ст)=бС1 si< ?)(f/+O И

6 6 е3т 2 (ст) = п-пст

р

„ . пв(пст) 1 , пст

+ Q sin — I — I + -ln--

66 1 2 l 2 ) 3 2

С + 1С1 cos ^)Р + O И)

Р-1 ~

2 l 2

d02 П п R„;„ пР (пст

= С1— Р sin du 2 2 l 2

-6+°М

dx2 1 „ п„ rcfifпст)^ 1 ,,/ r\

~¡b = 3ст + C1 ^ pcos п2в^п2ст^ +о(стР)

Отметим также, что вследствие того, что 9 -нечетная, т - четная по ст функции, значения

91 (0) = 0, ^^ (0) = 0. При подстановке этих соотно-

Сст

шений в (8) получим уравнение

пх

+

пст

П+с 2 С+е)со8?(т) + °И

- ЗС _3Х! (о)_1

Fr 2

2 ( 2

1 + ^ Г А (0)а_ПСТ +

2 2 ( dст 6

+ С БШ

2л 2

+ С» (а2в)Л

При приравнивании членов порядка 0(1) и 0(стр) получаются следующие уравнения

С _ т1 (о) - _

3

лр

3 1 л Рг' 1

решения задачи. Результаты представлены в логарифмическом масштабе. По оси ординат отложены десятичные логарифмы оценок относительных погрешностей _ ^ 5, 5-|ДРг/Рг| (точность полученных данных). Цифрой 0 отмечены оценки точности вычисленных данных, цифрами 1,2. ... - результаты первой, второй и т. д. фильтрации. Видно, что с помощью фильтрации удалось подавить несколько регулярных компонент (1) и на два порядка и более уменьшить погрешность искомого параметра Рг.

12

10

(1 + Р)С1^2 ,13

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

— - dz - dx + idy - е1^ - _е1ш(х^х И .

Для численного интегрирования применялась формула средних прямоугольников с последующей фильтрацией до достижения погрешности порядка 10-18.

Численное решение, обработка и анализ результатов. Численно задача решалась методом коллокаций. Уравнение (8) удовлетворялось в отдельных точках X - стт +10 (т - 0, п _ 1). Искомыми являлись значения Рг , С, 9т, т - 1, п _ 2. Полученная система п нелинейных уравнений решалась методом Ньютона с регулированием шага. Поиск решения прекращался, когда невязки по модулю не превышали 10-20. Для решения использовались числа с длиной мантиссы около 35 десятичных разрядов. Отличительной особенностью данного метода является использование при численном решении условия (8) вместо (4), что позволяет избежать численного интегрирования при формировании системы уравнений. Это существенно ускоряет решение системы. Решение задачи о солитоне Сток-са позволяет определить значения параметров волны, таких как число Фруда, амплитуда волны

а - ув _ 1 - Рг2/2 и т.д. В работе [13] была декларирована точность полученного значения числа Фруда порядка 13-ти разрядов. Позже в [1] был применен другой метод фильтрации, и точность доведена до 14-ти десятичных цифр. В [2] были применены другие методы вычисления и фильтрации, погрешность была декларирована на уровне 2 единиц 15-го разряда. Результаты последующих расчетов подтверждали оценки, полученные в предыдущих с помощью фильтрации (см. табл. 1).

На рис. 3а, 36 представлены результаты фильтрации данных, полученных этими методами

6 4 2

О

-1ВД

16 ■

12

I1« ■У**"" о

я •

* •

1 2 3 1ёи

а)

б

Рис. 3. Результаты фильтрации числа Фруда: а) согласно [13]; б) согласно [2]

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

Сравнение со значением, полученным в [2], показало отличие около _1,510-15. Тем самым, оценка, полученная ранее в [2], была подтверждена, и получено новое значение числа Фруда, равное Рг=1,2908904558633395±10-16. Кроме того, была модифицирована программная реализация используемого численного метода, что позволило значительно увеличить скорость вычисления.

Sherykhalina, AA. Sokolova // Mediterranean Journal of Soc. Sci. 2015. No. 2. P. 65-78.

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

4. Лаврентьев, МА. Методы теории функций комплексного переменного / М.А. Лаврентьев, Б.В. Ша-бат. - М.: Наука. 1973. 736 с.

5. Longuet-Higgins, M.S. On the mass, momentum, energy and circulation of the solitary wave / M.S. Longuet-Higgins, J.D. Fenton // Proc. R. Soc. Lond. 1974. V. 340. P. 471-493.

6. Fox, M.J.H. Ph.D. thesis. Cambridge Univ. 1977. 39 р.

7. Hunter, J.K Accurate computations for steep solitary waves / J.K. Hunter, J.-M. Vanden-Broek // J. Fluid. Mech. 1983. V. 136. P. 63-71.

8. Williams, J.M. Limiting gravity waves in water of finite depth. Phil. Trans. Roy. Soc. Lond. 1981. A302. N1466. P. 139-188.

9. Evans, W.A.B. An exact integral equation for solitary waves (with new numerical results for some internal properties). / W.A.B. Evans, M.J. Ford// Proc. R. Soc. London. 1996. A 452. P. 373-390.

10. Шерьхалина, Н.М. Разработка численных алгоритмов решения задач гидродинамики с особыми точками на свободной поверхности и экспериментальное исследование скорости их сходимости. Деп. в ВИНИТИ №2550-В95.

11. Maklakov, D.V. Almost-highest gravity waves on water of finite depth // Euro Journal of Applied Mathematics. 2002. 13. P. 67-93.

12. Житников, В.П. Расчет формы уединенных волн с помощью численно-аналитических методов / В.П. Житников, Н.М. Шерьхалина // Физика волновых процессов и радиотехнические системы. 1998. Т.1. № 2-3. С. 103-107.

13. 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 Journ. Vol. 10. No. 3. 2001. P. 309-314.

14. Гуревич, М.И. Теория струй идеальной жидкости. -М.: Наука, 1979. 536 с.

THE ERROR ESTIMATION AND ITS JUSTIFICATION BY THE FILTRATION OF NUMERICAL RESULTS OBTAINED WITH DIFFERENT NUMBERS OF NODES

© 2017 V.P. Zhitnikov, N.M. Sherykhalina, A.A. Sokolova

Ufa State Aviation Technical University

In the paper the obtained results error is estimated by the filtration of numerical data computed with different nodes numbers. In order to check the existence of possibility of programming errors and other additional errors the comparison of uncertainty intervals defined by different numerical methods is applied. The solution of the problem of Stokes's solitary wave is considered as an example. The boundary conditions are given in a differential form in this work. The solving of the problem by new method with greater accuracy allows establishing the fact of incorrectness of some estimates received earlier by the other authors and also say about reliability of the estimates received by the filtration methods.

Key words: numerical filtration, accuracy, calculation, Stokes's solitary wave

Vladimir Zhitnikov, Doctor of Physics and Mathematics,

Professor. E-mail: [email protected]

Nataliya Sherykhalina, Doctor of Technical Sciences.

E-mail: nsher @mail.ru

Alexandra Sokolova, Post-graduate Student.

E-mail: [email protected]

0 1 2 3 Ign

Рис. 4. Результаты фильтрации числа Фруда

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

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

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

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

2. Zhitnikov, V.P. Problem of Reliability Justification of Computation Error Estimates / V.P. Zhitnikov, N.M.