CC BY
17
Серия «Математика»
2015. Т. 11. С. 96—105
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 517.968
Применение численных методов для уравнений Вольтерра I рода, возникающих в обратной граничной задаче теплопроводности *
С. В. Солодуша
Институт систем энергетики им. Л. А. Мелентьева СО РАН
Аннотация. В прикладных задачах, связанных с исследованием нестационарных тепловых процессов, довольно часто возникает ситуация, когда невозможно осуществить прямые измерения требуемой физической величины и ее характеристики восстанавливаются по результатам косвенных измерений. При этом единственный путь отыскания требуемых значений связан с решением обратной задачи теплопроводности с исходными данными, известными только на части границы. Подобного рода задачи возникают не только при исследовании тепловых процессов, но и при исследовании процессов диффузии, изучении свойств материалов, связанных с тепловыми характеристиками. Статья посвящена численному решению уравнений Вольтерра I рода, полученных в результате применения интегральных преобразований Лапласа для решения уравнения теплопроводности. Работа состоит из введения и трех разделов. В первых двух разделах рассматривается специфика ядер из соответствующих интегральных уравнений. Указаны особенности их вычислений при выполнении операций машинной арифметики над вещественными числами с плавающей точкой. На тестах проиллюстрированы типовые случаи систематического накопления ошибок. В третьем разделе приведены результаты вычислительных алгоритмов, основанных на product integration method и квадратуре средних прямоугольников. Критериями выбора данных методов для решения интегральных уравнений Вольтерра I рода стали простота реализации вычислительной процедуры и возможность получения приближенного решения с погрешностью второго порядка по шагу сетки при точно заданных исходных данных. С целью проверки эффективности разностных методов приведены результаты тестовых расчетов.
Ключевые слова: обратная граничная задача теплопроводности, интегральные уравнения Вольтерра первого рода, численные методы.
* Работа выполнена при финансовой поддержке РФФИ, грант № 15-01-01425-а.
1. Введение
В статье [1] предложен способ поиска решения и(1,г) = ф(г), г > 0, обратной граничной задачи
иг = пхх, х е (0,1), г ^ 0, (1.1)
и (х, 0) = 0, и (0, г) = 0, их (0, г) = до (г) (1.2)
путем сведения (1.1), (1.2) к интегральному уравнению Вольтерра типа свертки:
г
Аф = ! к (г - эфв^ = у(г), 0 < 8 < г < т, (1.3)
о
те
к (г - 8) = £ (-1)9+1 д^-^-*, у (г) = ^ до (г). (1.4)
9=1
Задача (1.1), (1.2) играет важную роль в математических моделях из области нанотехнологий, метрологии, в изучении свойств материалов, связанных с тепловыми характеристиками [2; 3; 4; 5].
В работе [1], с учетом идей [6; 7], выполнена аппроксимация их(0,г) суммой N первых членов:
N
их(0,г) = 2п2^2 (-1)9+1 Я2 е-^^фШз,
9=1 о
где N — натуральное число. Тогда (1.3), (1.4) редуцируется к виду
j KN (г - 8)ф(8)й8 = у (г), (1.5)
о
N
KN(г - 8) = ^2 (-1)д+1 Я2е-п292(—). (1.6)
9=1
Обычно на практике вместо д0(г) известна (г): \д(г) - д0(г)\\с ^ 5, 5 > 0. Известно, что интегральные уравнения Вольтерра I рода относятся к классу условно-корректных задач, а процедура дискретизации обладает регуляризующим свойством с параметром регуляризации — шагом сетки, определенным образом связанным с уровнем возмущений исходных данных 5. Качественной теории и численным методам решения уравнений Вольтерра I рода посвящены многочисленные работы (см., например, [8; 9; 10; 11] и приведенную в них библиографию).
г
В данной статье рассматривается алгоритм приближенного решения (1.5), (1.6) при точно заданной правой части. Заметим, что при решении (1.5), (1.6) возникают три вида погрешностей, связанных, во-первых, с аппроксимацией исходной задачи (1.3), (1.4), во-вторых, с погрешностью метода и, наконец, с вычислительными ошибками при выполнении операций машинной арифметики над вещественными числами с плавающей точкой. Наиболее интересен для исследования первый из указанных случаев. Однако для перехода к проблеме оценивания параметра N в (1.6) следует разработать алгоритм вычисления Кмтах, учитывающий особенности машинной арифметики и обеспечивающий желаемое (заданное) число достоверных позиций в мантиссе.
Цель данной статьи — рассмотреть применение численных методов для уравнений вида (1.5), (1.6) с учетом механизмов возникновения погрешностей при компьютерных расчетах. Работа продолжает исследование, начатое в [12].
2. Свойства ядер (1.6)
те
П2 ^ Км(в)(1в =
Для понимания специфики (1.5), (1.6) рассмотрим свойства ядер Км при фиксированных значениях N. Очевидно, что для любых N имеет место
те
1, если N нечетно, 0, если N четно,
о
так как
те
(_1)м+1 п2N2 [ е-п2м2Чз = ( 1, если N нечетно,
у ^ _1, если N четно.
о
Пусть Л* есть решение уравнения Км (Л) = 0, где Л = £ _ в. Тогда, по формуле Ньютона-Лейбница:
/Л 1 (_1)М+1 м
КмШв = 2 + _^ + С1, С1 = £ (-1Г д2е-п^, 0 4=1
где С1 с ростом N стремится к — 1.
Значения Км на диагонали £ = в при N = 10,21, а также Л*, полученные с одинарной точностью, приведены в таблице 1.
Вычислительный эксперимент в [12] показал, что с ростом числа слагаемых в (1.6) значения Л* убывают (при этом имеется монотонность только отдельно по четным и нечетным N). Из таблицы видно, что Км(0) = 0. Как известно, это условие является необходимым для
корректности (1.5), (1.6) на паре (С,С^ Т]), где у(0) = 0, у'(1) е С[0,Т]. Далее используем Л* для ограничения величины шага сетки Н сверху, чтобы значение сеточной функции КN в первом узле было отличным от нуля.
Таблица 1
Числовые характеристики ядер КN.
N Л* Kn (0) N Л* Kn (0) N Л* Kn (0)
10 0.01378 -542.828 14 0.01019 -1036.308 18 0.00809 -1687.702
11 0.01221 651.394 15 0.00789 1184.353 19 0.00516 1875.225
12 0.01173 -769.829 16 0.00913 -1342.266 20 0.00735 -2072.617
13 0.01022 898.134 17 0.00631 1510.049 21 0.00429 2279.879
3. Типовые случаи возникновения погрешностей при компьютерных расчетах
Рассмотрим известные случаи систематического накопления ошибок (см., например, [13]), возникающие при вычислении значений ядра (1.6). Используем систему компьютерной алгебры Maple 10. Следуя [14], включим в общепринятое представление вещественного числа параметр f, равный числу достоверных позиций в мантиссе (начиная слева). Считаем, что x = s ■ M ■ 10-m+p задается набором (s,M,p,f), где s £ {-1, 0, +1} — знак числа, M £ {10т-1, 10т-1 + 1,..., 10т - 1}U {0} — мантисса числа, m — число позиций мантиссы, p — порядок числа.
Проиллюстрируем специфику расчетов при сложении чисел разных порядков в (1.6) на примере.
Пример 1. Пусть N = 50, Ао = 10-3, m > 8. Выберем
34 50
X! = £ (-1)q+1 д2в-п212х°, Х2 =^2 (-1)q+1 q2e-n2q2x° q=11 q=35
и найдем xs = x1 + x2.
В таблице 2 даны параметры (1, M1, 2, f1), (1, M2, -2, f2) и (1, Ms, 2, fs), которые определяют xi, x2 и xs соответственно.
Следующий пример иллюстрирует ситуацию, возникающую при вычислении разности между числами в (1.6), у которых совпадают порядки и несколько старших позиций мантиссы.
Пример 2. Пусть N = 50, А0 = 10-3, m > 8. Введем
10 50
x3 = £ (-1)q+1 g2e-n2q2x°, x4 =J2 (-1)q+1 q2e-n2q2Xo. q=1 q=11
Определим хд = \х4\ — \х3\.
Таблица 2
Значения М и / для х1, х2 и ху.
т М1 /1 М2 /2 Му /у
8 18652239 6 44981421 6 18656737 6
9 186522441 8 449814458 7 186567422 8
10 1865224455 9 4498144699 8 1865674269 8
11 18652244592 11 44981446726 8 18656742737 11
12 186522445926 11 449814466957 10 186567427373 11
13 1865224459248 12 4498144669376 10 1865674273715 12
14 18652244592468 13 44981446694089 11 18656742737137 13
15 186522445924647 15 449814466941161 13 186567427371341 14
Ниже даны параметры (-1, М3, 2, /3), (1, М4, 2, /4) и (-1, Мд, 2, /д), которые задают значения ж3, х4 и хд.
Таблица 3
Значения М и / для хз, х4 и хд.
т Мз /з М4 /4 Мд /д
8 18656743 7 18656737 6 00000006 0
9 186567428 8 186567424 8 000000004 0
10 1865674274 9 1865674268 8 0000000006 0
11 18656742750 11 18656742736 10 00000000014 1
12 186567427505 12 186567427372 11 000000000133 3
13 1865674275054 12 1865674273715 12 0000000001339 4
14 18656742750534 14 18656742737138 13 00000000013396 4
15 186567427505342 14 186567427371340 14 000000000134002 2
Очевидно, что обнуление нескольких старших позиций приводит к возникновению числа с меньшим количеством значащих цифр в мантиссе.
Покажем на тестовом примере важность согласования количества цифр т с числом слагаемых N в (1.6).
Пример 3. Пусть в (1.6) N = 10,100. Предположим, что гарантированное число позиций мантиссы /5 = 8 (/б = 16). Определим соответствующие значения т1 (т2), обеспечивающие достоверные значения Км (А,), г = 1,2, при А1 = 10"2, А2 = 10"3.
Результаты расчетов проиллюстрированы на рисунке 1. Сплошной линией изображен график т1, которому отвечает желаемое значение /5.
Пунктирной линией дан график которому соответствует значение /б. С правой стороны отмечены Х^, при которых проводились расчеты К^- При подсчете КN(Х2) было получено (см. рис. 1), что с увеличением N требуемая точность обеспечивается за счет существенного роста т. Дальнейшее вычисление KN(Х2) при N > 100 показало постоянную среднюю скорость роста значений т.
Учитывая перечисленные механизмы возникновения вычислительной ошибки при компьютерных расчетах, перейдем теперь к задаче приближенного решения (1.5), (1.6).
4. Результаты приближенного решения (1.5), (1.6)
Введем на [0, T] равномерные сетки узлов ti = ih, ti-1 = (i — 2)h, i =
1, n, nh = T и, аппроксимируя интеграл в (1.5) с помощью квадратуры средних прямоугольников и product integration method [15], запишем соответствующие сеточные аналоги (1.5), (1.6)
N i
72У" e-n2q2(j-2)h jh.
h£ ( —1)q+1 q2£ e-n q U-1 )h 1 = y q=1 j=1 2
n i jh
E(—1)q+1 q2E j| / e- ih-s')ds = yh.
q=1 j=1 (j-1)h
и
Обозначим их решение через <ph и <fih соответственно. Проведем вычислительный эксперимент, обеспечивающий заданное число /5 = 8 достоверных цифр в мантиссе.
Пример 4. В качестве модельной функции выберем ф (t) из [16]:
_t
-г / ч 1 — e с
Ф (t) =-1 — t,
1 — e-1
где е = 10-1. Зададим N = 2, 5.
В табл. 4 приведены результаты численных расчетов на отрезке интегрирования [0,1]. Здесь приняты следующие обозначения: ||eh||cfe — максимум модуля отклонений в узловых точках точного решения от приближенного, полученного с помощью квадратуры средних прямоугольников; ||eh||cfe — максимум модуля отклонений в узловых точках точного решения от приближенного, полученного с помощью product integration method.
Таблица 4
Погрешности сеточного решения.
h р->< QJ; to m> OS; l to II-||N=3 1 1 £ 1 1 Ch CN yll CO
1/128 0.020527 0.001985 0.015342 0.004622
1/256 0.005001 0.000499 0.003815 0.001171
1/512 0.001242 0.000125 0.000952 0.000294
1/1024 0.000310 0.000031 0.000238 0.000074
h 1 1*1 1 N= ll-l |N=4 1 1 e 1 1 Ch ||-||N=5 1 1 e 1 1 Ch ||-||N=6 1 1 £ 1 1 Ch
1/128 0.366901 0.007947 0.105599 0.011669
1/256 0.006501 0.002056 0.025682 0.003130
1/512 0.015090 0.000519 0.006377 0.000797
1/1024 0.003707 0.000129 0.001591 0.000200
Из таблицы следует, что оба разностных метода имеют порядок сходимости O(h2), хотя при сравнении ||eh||cfe и ||eh||cfe видно, что применение метода product integration предпочтительней.
Замечание 1. Вычислительный эксперимент, проведенный для данного примера при е = 10-3, показал сходимость для h < 10-3. В связи с этим в дальнейшем планируется привлечение численных методов более высокого порядка, в частности, методов Рунге-Кутта третьего и четвертого порядков.
5. Заключение
Подводя итог, отметим, что в данной статье исследовано приближенное решение интегрального уравнения Вольтерра I рода типа свертки, возникающего в обратной граничной задаче теплопроводности, с помощью разностных методов второго порядка. Приведены результаты расчетов с использованием системы Maple 10. Вычислительный эксперимент проводился с учетом механизмов возникновения погрешностей при компьютерных расчетах. На тестовых примерах показаны типовые случаи систематического накопления ошибок.
Автор выражает благодарность А.С. Апарцину и И.В. Мокрому за оказанную помощь при обсуждении результатов.
Список литературы
1. Yaparova N. M. Numerical methods for solving a boundary value inverse heat conduction problem / N. M. Yaparova // Inverse Problems in Science and Engineering. - 2014. - Vol. 22, N. 5. - P. 832-847.
2. Идентификация математических моделей теплопереноса в разлагающихся материалах / О. М. Алифанов, С. А. Будник, А. В. Ненарокомов, А. В. Нетелев // Тепловые процессы в технике. - 2011. - № 8. - С. 338-347.
3. Модель роста нанокристаллов в аморфном сплаве / П. А.Гамов, А. Д. Дрозин, М. В. Дудоров, В. Е. Рощин // Металлы. - 2012. - № 6. - С. 101-106.
4. Короткий А. И. Реконструкция граничных режимов в обратной задаче тепловой конвекции высоковязкой жидкости / А. И. Короткий, Д. А. Ковтунов // Тр. Ин-та математики и механики УрО РАН. - 2006. - Т. 12, № 2. - С. 88-97.
5. Балабанов П.В. Математическое моделирование теплопереноса в процессе хе-мосорбции / П. В. Балабанов, С. В. Пономарев, А. В. Трофимов // Вестн. ТГТУ. - 2008. - Т. 14, № 2. - С. 334-341.
6. Beilina L. Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems / L. Beilina, M.V. Klibanov. - N. Y. : Springer, 2012.
7. Kabanikhin S. I. Inverse and Ill-Posed Problems. Theory and Applications / S. I. Kabanikhin. - Germany : De Gruyter, 2011.
8. Brunner H. The Numerical Solution of Volterra Equations / H. Brunner, P. J. van der Houwen. - North-Holland, Amsterdam, 1986.
9. Brunner H. Collocation methods for Volterra integral and related functional differential equations / H. Brunner. - N. Y. : Cambridge Univ. Press, 2004.
10. Верлань А. Ф. Интегральные уравнения: методы, алгоритмы, программы / А. Ф. Верлань, В. С. Сизиков. - Киев : Наук. думка, 1986.
11. Апарцин А. С. Неклассические уравнения Вольтерра I рода: теория и численные методы / А. С. Апарцин. - Новосибирск : Наука, 1999.
12. Solodusha S. V. Numerical solution of the Volterra equations of the first kind that appear in an inverse boundary-value problem of heat conduction / S. V. Solodusha, N. M. Yaparova // to appear in Numerical Analysis and Applications, http://arxiv.org/abs/1407.1678 (date of access: 01.11.2014).
13. Калиткин Н. Н. Численные методы / Н. Н. Калиткин. - М. : Наука, 1978.
14. Мокрый И. В. Основные механизмы возникновения вычислительной ошибки при компьютерных расчетах / И. В. Мокрый, О. В. Хамисов, А. С. Цапах
// Материалы IV Всерос. конф. «Проблемы оптимизации и экономические приложения». - Омск : Наследие, 2009. - С. 185.
15. Linz P. Product integration method for Volterra integral equations of the first kind / P. Linz // BIT. - 1971. - Vol. 11. - P. 413-421.
16. Geng F. Z. Analytical Approximation to Solutions of Singularly Perturbed Boundary Value Problems / F. Z. Geng, M. G. Cui // Bulletin of the Malaysian Mathematical Sciences Society. -2010. - Vol. 33, N 2. - P. 22-232.
Солодуша Светлана Витальевна, кандидат физико-математических наук, доцент, Институт систем энергетики им. Л. А. Мелентьева СО РАН, 664033, Иркутск, ул. Лермонтова, 130, тел.: (3952)500646 (e-mail: [email protected])
S. V. Solodusha
Application of Numerical Methods for the Volterra Equations of the First Kind that Appear in an Inverse Boundary-Value Problem of Heat Conduction
Abstract.In applied problems related to the study of non-stationary thermal processes, often arise a situation where it is impossible to carry out direct measurements of the desired physical quantity and its characteristics are restored on the results of indirect measurements. In this case the only way to finding the required values is associated with the solution of the inverse heat conduction problem with the initial data, known only to the part of the boundary. Such problems appear not only in the study of thermal processes, but also in the study of diffusion processes and studying the properties of materials related to the thermal characteristics. This article is devoted to the approximate solution of the Volterra equations of the first kind received as a result of the integral Laplace transform to solve the heat equation. The work consists of an introduction and three sections. In the first two sections the specificity of Volterra kernels of the corresponding integral equations and peculiarity of computing kernels over the machine arithmetic operations on real numbers with floating point are considered. In tests typically systematic accumulation of errors are illustrated. The third section presents the results of numerical algorithms based on product integration method and middle rectangles quadrature. The conditions under which used algorithms are stable and converge to the exact solution in the case of fixed digit grid in the computer representation of numbers are allocated. Series of test calculations are carried out in order to test the efficacy of difference methods.
Keywords: inverse boundary-value problem of heat conduction, Volterra integral equations of the first kind, numerical methods.
References
1. Yaparova N.M. Numerical methods for solving a boundary value inverse heat conduction problem. Inverse Problems in Science and Engineering, 2014, vol. 22, no 5, pp. 832-847.
2. Alifanov O.M., Budnik S.A., Nenarokomov A.V., Netelev A.V. Identification of Mathematical Models of Heat Transfer in Decomposing Materials [Identifikatsiya
matematicheskikh modeley teploperenosa v razlagayushchikhsya materialakh]. Teplovye protsessy v tekhnike [Thermal processes in engineering], 2011, no 8, pp. 338-347.
3. Gamov P.A., Drozin A.D., Dudorov M.V., Roshchin V.E. The Growth Model of Nanocrystals in the Amorphous Alloy [Model' rosta nanokristallov v amorfnom splave]. Metally [Metals], 2012, no 6, pp. 101-106.
4. Korotkii A.I., Kovtunov D.A. Reconstruction of the boundary conditions in the inverse problem of thermal convection of a viscous fluid [Rekonstruktsiya granichnykh rezhimov v obratnoy zadache teplovoy konvektsii vysokovyazkoy zhidkosti]. Proceedings of the Steklov Institute of Mathematics, 2006, vol. 255, no 2, pp. S81-S92.
5. Balabanov P.V., Ponomarev S.V., Trofimov A.V. Mathematical Modeling of Heat Transfer During of Chemosorption [Matematicheskoe modelirovanie teploperenosa v protsesse khemocorbtsii]. Vestnik TGTU [Bulletin of the Tambov State Technical University], 2008, vol. 14, no 2, pp. 334-341.
6. Beilina L., Klibanov M.V. Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems. New York, Springer, 2012.
7. Kabanikhin S.I. Inverse and Ill-Posed Problems. Theory and Applications. Germany, De Gruyter, 2011.
8. Brunner H., van der Houwen P.J. The Numerical Solution of Volterra Equations. North-Holland, Amsterdam, 1986.
9. Brunner H. Collocation methods for Volterra integral and related functional differential equations. New York, Cambridge Univ. Press, 2004.
10. Verlan' A.F., Sizikov V.S. Integral equations: methods, algorithms, programs [Integralnye uravneniya: metody, algoritmy, programmy]. Kiev, Nauk. dumka, 1986.
11. Apartsyn A.S. Nonclassical linear Volterra equations of the first kind. Boston, VSP Utrecht, 2003.
12. Solodusha S.V., Yaparova N.M. Numerical solution of the Volterra equations of the first kind that appear in an inverse boundary-value problem of heat conduction. To appear in Numerical Analysis and Applications, available at: http://arxiv.org/abs/1407.1678 (accessed 1 November 2014).
13. Kalitkin N.N. Numerical methods [Chislennye metody]. M., Nauka, 1978.
14. Mokry I.V., Khamisov O.V., Tsapakh A.S. The Basic Mechanisms of the Emergence of Computational Errors in Computer Calculations [Osnovnye mekhanizmy vosniknoveniya vychislitelnoy oshibki pri kompyuternykh raschetakh]. Materialy IV Vseross. konf. "Problemy optimizatsii i economicheskie prilozheniya" [Proc. IVth All-Russian Conference "Problems of Optimization and Economic Applications"], Omsk, Nasledie, 2009, pp. 185.
15. Linz P. Product integration method for Volterra integral equations of the first kind. BIT, 1971, vol. 11, pp. 413-421.
16. Geng F.Z., Cui M.G. Analytical Approximation to Solutions of Singularly Perturbed Boundary Value Problems. Bulletin of the Malaysian Mathematical Sciences Society, 2010, vol. 33, no 2, pp. 221-232.
Solodusha Svetlana Vital'evna, Associate Professor, Melentiev Energy Systems Institute SB RAS, 130, Lermontov st., Irkutsk, 664033, tel.: (3952)500646 (e-mail: [email protected])
