CC BY
23Математические заметки СВФУ Январь—март, 2019. Том 26, № 1
УДК 519.63
ИТЕРАЦИОННАЯ ИДЕНТИФИКАЦИЯ СТАЦИОНАРНОЙ ПРАВОЙ ЧАСТИ ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ Л. Д. Су, В. И. Васильев
Аннотация. Для одномерного параболического уравнения рассмотрена задача определения правой части, зависящей только от пространственных переменных. Для численного решения поставленной обратной начально-краевой задачи используется метод сопряженных градиентов в сочетании с методом конечных разностей с неявной аппроксимацией по времени с весовым множителем а £ [0, 1]. Возможности предложенного вычислительного алгоритма подтверждены результатами вычислительного эксперимента для модельных задач с квазиреальными решениями, включая задачи с условиями переопределения, имеющими случайные ошибки.
Б01: 10.25587/8УРи.2019.101.27249
Ключевые слова: параболическое уравнение, обратная задача, идентификация источника, разностная схема, метод сопряженных градиентов.
1. Введение
Решение многих прикладных актуальных задач науки и техники сводится к обратным задачам. Вопросы их теоретического исследования и разработка эффективных численных методов решения различных обратных задач математической физики обобщены в монографиях А. А. Самарского, П. Н. Вабищевича [1], С. И. Кабанихина [2] и Лионса, Мадженеса [3], В. Исакова [4], К. Фогеля [5].
С точки зрения приложений большой интерес представляют обратные задачи идентификации источника в уравнении теплопроводности. Различным аспектам исследования вопросов корректности постановки обратных задач, в том числе и задач идентификации пространственно-распределенного источника для параболических уравнений, посвящены работы [6-9].
Йохансон и Лесник в [10,11] предложили итерационный алгоритм устойчивого итерационного восстановления источника с использованием метода граничных элементов, в работе [12] — метода конечных элементов, в [13] метода конечных разностей. В этих методах на каждой итерации решается корректная прямая задача. П. Н. Вабищевич в работе [14] развил, обосновал и численно реализовал предложенный им ранее в монографии [1] итерационный метод и
Работа выполнена при поддержке гранта Российского фонда фундаментальных исследований (проект № 17-01-00732 А).
© 2019 Су Л. Д., Васильев В. И.
построил второй итерационный процесс типа Пикара. Оба итерационных метода не содержат параметров регуляризации. В работе [15] предложен итерационный метод вариационного типа, основанный на использовании множителей Лагранжа для идентификации оптимальных значений параметров итерационного процесса. В работе А. Л. Карчевского [16] предложен метод сопряженного оператора решения обратных задач подобного типа.
Следует отметить, что в [17] для численного решения ретроспективной обратной задачи теплопроводности был предложен итерационный метод минимальных невязок, а в [18] — метод сопряженных градиентов.
В данной работе для численного решения задачи идентификации источника, зависящего от пространственных переменных в параболическом уравнении, предлагается использовать итерационный метод сопряженных градиентов. Метод достаточно прост и эффективен. Приведены примеры решения модельных тестовых задач. Представлены результаты расчетов при наличии шума.
2. Постановка задачи
В области От = О х (0, Т], О = [0,1] рассмотрим обратную задачу для линейного параболического уравнения с переменными коэффициентами. Ищем пару функций и(х,г), /(х) из условий: искомая функция и(х,г) в области От является решением параболического уравнения второго порядка
и удовлетворяет граничным и начальному условиям:
и(ж, г) = о, х е до, о < г < т, (2.2)
и(х,0) =и0(х), же П. (2.3)
Пусть для идентификации неизвестной функции /(х) задано условие переопределения
и(х,Т) = <р(х), х£П. (2.4)
Если коэффициенты уравнения (2.1) достаточно гладкие функции, удовлетворяющие условиям 0 < С1 < к(х,г) < с2 < ж, функции д(ж,г), «о(х) достаточно гладкие и ограниченные, кроме того, начальное условие обращается в нуль на границе области О, то обратная задача корректна при соответствующих условиях существования и однозначной разрешимости, которые можно найти в [4,6,8]. В данной работе мы основное внимание уделяем проблемам численного решения рассматриваемых обратных задач, не останавливаясь на теоретических вопросах сходимости приближенного решения к точному.
3. Разностный аналог задачи
Численное решение параболической задачи (2.1)—(2.4) проведем с помощью конечно-разностного метода [19]. Для этого в области О введем равномерную сетку с постоянным шагом к и обозначим через ш множество внутренних узлов:
ш = = гк, г = 1,. .., п — 1}.
Введем гильбертово пространство сеточных функций у,у € Н = Ь2(ш), в котором скалярное произведение и норма определены следующим образом:
п-1
(у, = У1Щ,г> =
г=1
В предположении достаточной гладкости коэффициента к(х, /) для всех внутренних узлов х € ш сеточный аналог оператора А(£) запишем в виде
Л(г)у = -(а(х,г)ух)х, х € ш. (3.1)
Здесь в соответствии с интегро-интерполяционным методом построения дискретного аналога оператора (3.1) можно пользоваться формулами вычисления коэффициентов:
а(х, /) = к(х + 0.5/г, /), а(х, /) = — (к(х, /) + к(х + /г., £)),
а(ж, = ^(1//г(ж, + 1/к(х + /г,
В множестве сеточных функций, обращающихся в нуль на множестве граничных узлов дш, построенный оператор А(£) является самосопряженным, положительно определенным и ограниченным:
8п 4го
Обратной задаче (2.1)-(2.4) поставим в соответствие задачу Коши для дифференциально-операторного уравнения
Г ^ + = /(ж)д(ж,/), ж е о < г < т,
{, у(х, 0) = щ(х), х £ ш,
при заданном дополнительном условии переопределения
у(х, Т) = <р(х), х£ш. (3-3)
Обозначим через ут разностное решение на момент времени = тт, где т > 0 — шаг по времени, причем Мт = Т. Дифференциально-разностной задаче (3.2), (3.3) поставим в соответствие разностную схему с весовым множителем а € [0,1]:
ут+1 _ ут
У--—+ Ат+а(аут+1 + (1 - а)ут) = /(х)д(х, 1т + *т), 4)
х € ш, т = 0,1,. .., М - 1,
уо = щ(х), х £ ш. (3-5)
Здесь при решении обратной задачи идентификации дискретного аналога пространственно-распределенного источника /(ж), х € и, используем сеточный аналог дополнительного условия (3.3)
уК(х) = (р(х), х £ и. (3-6)
4. Алгоритм решения обратной задачи
При приближенном решении разностной обратной задачи (3.4)—(3.6) для уточнения начального условия будем пользоваться трехслойным итерационным методом сопряженных градиентов. При этом на каждой итерации решаются корректные прямые задачи с помощью стандартных двуслойных разностных схем [19].
Придадим этой задаче соответствующую операторную формулировку. Из (3.4), (3.5) для заданного у0 = ио(х) на конечный момент времени получим
м м
ум = (р = я/щ + «{/ = Бт, ?М=^тсёт, (4.1)
т=1 т=1
где Бт — оператор перехода с одного т-го временного слоя на следующий (т + 1)-й временной слой:
Бт = (1 + агЛт+а)-1(1 - (1 - <г)тЛт+а), (4.2)
а оператор определяется по формуле
т
Vт = П д(х, + ^т )(1 + ТаЛт+а )-1. (4.3)
3=1
Таким образом, приближенное решение обратной задачи (2.1)—(2.4) с помощью разностной схемы (3.4)—(3.6) приводит к решению следующего сеточного операторного уравнения:
8§/(х) = ф(х)-£/щ, хеи. (4.4)
В силу самосопряженности оператора Л самосопряженным является оператор с£п и, следовательно, оператор 8$ в уравнении (4.4). Однозначная разрешимость операторного уравнения (4.4) следует из положительности оператора 83.
Для уменьшения количества итераций целесообразно ориентироваться на применение итерационных методов вариационного типа. При использовании наиболее быстро сходящегося итерационного метода сопряженных градиентов имеем следующий итерационный процесс [20].
1. Полагаем к = 0.
2. Следуя П. Н. Вабищевичу [21], начальное приближение искомой правой задаем с помощью формулы
/о(ж) = ~(а(х, Т)1рх)х/д{х, Т), х£ш.
Определяем начальное приближение решение дискретного аналога прямой задачи для параболического уравнения с заданным начальным условием и начальным приближением правой части
Ув=щ(х), хеш,
m+1
---+ + (1=fo{x)g{x,tm + (TT), (45)
x £ ш, m = 0,1,..., M - 1,
2. Вычисляем начальную невязку го(х) = уq (ж) — <р(х), х £ си, и задаем начальное приближение вспомогательной сеточной функции po(x) = ro(x), x £ ш.
3. Запускаем счетчик итераций k = k + 1.
4. Решаем дискретный аналог прямой задачи для вспомогательного параболического уравнения:
= 0, х £ U,
Zm+1 _ zm
^-^ + Arn+a{azlL+1 + {1 - = pkg{x,tm +(тт), (4 g)
x £ ш, m = 0,1, ..., M - 1,
5. Определяем новое значение первого итерационного параметра
(г к, г к) , ,
ак = т-г. 4.7
(Zk,Pk)
6. Вычисляем очередные приближения искомой правой части и невязки по формулам
fk+1 = fk + akPk, rk+i = rk - akZk, x £ ш. (4.8)
7. Находим очередной итерационный параметр и значение вспомогательного вектора
Рк = ^Гк+1,Гк+1\ Рк+1 = Гк + РкРк, х £ ui. (4.9)
(rk ,rk)
8. Итерационный процесс продолжается до тех пор, пока не выполнится критерий остановки итераций \\rk||/|М1 < е, в противном случае продолжаем процесс, возвращаясь к п. 4.
5. Численные расчеты
В этом разделе приводятся результаты вычислительного эксперимента, подтвердившие эффективность предложенного итерационного метода численного решения обратной задачи идентификации стационарного источника параболического уравнения. Введем относительные ошибки в сеточных пространствах C, 12 следующим образом:
тах|Д(Ж)-/(Ж)| \\fk(x)-f(x) ||
кс —-;—| г, -, щ2 —
max |f (x)| ' 12 \\f (x)
где /к (х) — к-е приближение функции источника, /(х) — точное (квазиреальное) распределение источника.
Поскольку имеем дело с обратной задачей, в которой задается дополнительное условие, а оно, как правило, известно с некоторой случайной ошибкой
(«шумом»). Шум зададим с помощью генератора случайных чисел, равномерно распределенных в сегменте [-0.5, 0.5]:
tps (x) = ^(x) + (5.1)
Пример 1. Пусть искомый стационарный источник задан достаточно гладкой функцией, имеющей вид нормального распределения Гаусса:
f (ж)= e-(x-0-51)2, x G [0,1]. (5.2)
Расчеты проводились при следующих значениях исходных данных: 1 = 20; T =1; N = M = 1000; h = 1/N; т = T/M;
k(x,t) = 1, g(x,t) = 1, (x,t) G [0,1] x [0,T]; u0(x) = 0, x G [0,1].
Условие переопределения получаем из решения прямой задачи с известной правой частью (5.2). Решение прямой задачи осуществляем с помощью разностной схемы повышенного порядка точности с весовым множителем <г = 0.5 — h2/(12T) с порядком аппроксимации 0(т2 + h4) на достаточно подробной пространственно-временной сетке, и ее, следуя П. Н. Вабищевичу, назовем «квазиреальным» значением прямой задачи решением в финальный момент времени.
На рис. 1 (слева) показан график решения прямой задачи в финальный момент времени, справа — график заданной правой части f (ж), и здесь точками представлены значения идентифицированной правой части. Реализация предложенного итерационного метода осуществлялась на достаточно грубой сетке с N = M = 200, h = 1/N, т = T/M с условием выхода из итерационного цикла
1Ы1/1М1 <£, £ = 10-7.
На рис. 2 представлены графики относительной ошибки рассматриваемого итерационного метода решения поставленной задачи на разных пространственно-временных сетках с N = 200: слева — в C, справа — в L2. Видно, что итерационный процесс сходится очень быстро и с увеличением числа разбиений точность разностной схемы повышается.
Результаты счета для обратной задачи заданным «шумом» с разными значениями параметра S = 0.01, 0.03 представлены на рис. 3. Видно, что с ростом уровня «шума» погрешность идентификации правой части растет.
Пример 2. Будем искать стационарную правую часть в виде разрывной функции, заданной по формуле
{1, ж G (0.41, 0.61), 0.5, ж = 0.41&0.61, 0, в противном случае.
График решения прямой задачи в финальный момент времени представлен на рис. 4 (слева), а справа — график заданной правой части f(x) на очень
Рис. 2. Графики относительных ошибок на разных сетках Не (слева) и Нь2 (справа) на каждой итерации гладкой правой части
Рис. 3. Графики вычисленной гладкой правой части при разных уровнях «шума»
подробной сетке с N = М = 1000. Точками представлены значения идентифицированной правой части. Здесь реализация предложенного итерационного метода также осуществлялась на достаточно грубой сетке с N = М = 200.
На рис. 5 представлены графики относительной ошибки рассматриваемого итерационного метода решения поставленной задачи на разных пространственно временных сетках с N = 200: слева — в С, справа — в Видно, что
0.6 ^ 0.4 0.2
- Ехай * » Митепса!
Рис. 4. Графики функций <р(ж) (слева) и разрывной правой части /(ж) (справа).
ю1
в? ю-3
• • М=100 ■ М=200 • о М=300 ■ ■ М=400
оГ ю
• • м=юо
• • Л/= 200
М = 300
■ • Л/=400
• • • • •
01234567 к
01234567 к
Рис. 5. Графики относительных ошибок на разных сетках Не (слева) и Нь2 (справа) на каждой итерации идентификации разрывной правой части
1.0
1.0
0.8
0.8
0.6
0.4
0.2
0.0
0.0
и
5
-1
-2
-5
Рис. 6. Графики вычисленной разрывной правой части при разных уровнях «шума»
итерационный процесс сходится очень быстро и с увеличением числа разбиений точность разностной схемы повышается.
Результаты счета для обратной задачи заданным «шумом» с разными значениями параметра 6 = 0.01, 0.05 представлены на рис. 6. Видно, что с ростом уровня «шума» погрешность идентификации правой части растет.
6. Заключение
В работе предложен итерационный метод решения обратной задачи идентификации стационарного источника тепла. Приведенные результаты вычислительного эксперимента подтвердили высокую эффективность итерационного метода сопряженных градиентов для численного решения обратной задачи идентификации правой части параболического уравнения, зависящего только от пространственных переменных.
Благодарности
Авторы выражают искреннюю благодарность профессору П. Н. Вабищеви-чу за конструктивные замечания и плодотворные обсуждения.
ЛИТЕРАТУРА
1. Samarskii A. A., Vabishchevich P. N. Numerical methods for solving inverse problems of mathematical physics. Berlin: De Gruyter, 2007.
2. Kabanikhin S. I. Inverse and ill-posed problems. Theory and applications. Berlin: De Gruyter, 2011.
3. Lions J.-L., Magenes E. Non-homogeneous boundary value problems and applications. Berlin: Springer-Verl., 1972. V. I.
4. Isakov V. Inverse source problems. Providence, RI: Amer. Math. Soc., 1990. (Math. Surv. Monogr.; V. 34).
5. Vogel C. R. Computational methods for inverse problems. Society for Industrial and Applied Mathematics, 2002.
6. Cannon J. R. Determination of an unknown heat source from overspecified boundary data // SIAM J. Numer. Anal. 1968. V. 5, N 2. P. 275-286.
7. Isakov V. Inverse parabolic problems with the final overdetermination // Commun. Pure Appl. Math. 1991. V. 44, N 2. P. 185-209.
8. SoloVev V. V. Solvability of the inverse problems of finding a source, using overdetermination on the upper base for a parabolic equation // Diff. Equ. 1990. V. 25. P. 1114-1119.
9. Kostin A. B. The inverse problem of recovering the source in a parabolic equation under a condition of nonlocal observation // Sbornik Mathematics. 2013. V. 204. P. 1391-1434.
10. Johansson T., Lesnic D. Determination of a spacewise dependent heat source //J. Comput. Appl. Math. 2007. V. 209, N 1. P. 66-80.
11. Johansson B. T., Lesnic D. A variational method for identifying a spacewise-dependent heat source // IMA J. Appl. Math. 2007. V. 72, N 6. P. 748-760.
12. Dhaeyer S., Johansson B. T., Slodicka M. Reconstruction of a spacewise- dependent heat source in a time-dependent heat diffusion process // IMA J. Appl. Math. 2014. V. 79, N 1. P. 33-53.
13. Erdem A., Lesnic D., Hasanov A. Identification of a spacewise dependent heat source // Appl. Math. Model. 2013. V. 37, N 24. P. 10231-10244.
14 Vabishchevich P. N. Additive operator-difference schemes. Berlin; Boston: Walter de Gruyter GmbH, 2014.
15. Liu J., Tang J. Variational iteration method for solving an inverse parabolic equation // Phys. Let. A. 2008. V. 372, N 20. P. 3569-3572.
16. Karchevsky A. L. Reformulation of an inverse problem statement that reduces computational costs // Euras. J. Math. Comp. Appl. 2013. V. 1, N 2. P. 4-20.
17. Samarskii A. A., Vabishchevich P. N., Vasil'ev V. I. Iterative solution of a retrospective inverse problem of heat conduction // Matematicheskoe Modelirovanie. 1997. V. 9, N 5. P. 119-127.
18. Vasil ev V. I., Kardashevsky A. M. Iterative solution of the retrospective inverse problem for a parabolic equation using the conjugate gradient method // Lect. Notes Comput. Sci. 2-7. V. 1087. P. 666-673.
19. Samarskii A. A. The theory of difference schemes. New York: Marcel Dekker, 2001.
20. Saad U. Iterative methods for sparse linear systems. 2-nd Edition. SIAM, 2003.
21. Vabishchevich P. N. Iterative computational identification of a spacewise dependent the source in a parabolic equations // Inverse Probl. Sci. Engineering. 2016. P. 1—23.
Поступила в редакцию 10 февраля 2019 г. После доработки 26 февраля 2019 г. Принята к печати 1 марта 2019 г.
Су Лин Де, Васильев Василий Иванович
Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики, ул. Кулаковского, 42, Якутск 677000 sulingde@gmail•com, [email protected]
Математические заметки СВФУ Январь—март, 2019. Том 26, № 1
UDC 519.63
ITERATIVE IDENTIFICATION OF THE SPACEWISE-DEPENDENT RIGHT-HAND SIDE IN A PARABOLIC EQUATION L. D. Su and V. I. Vasil'ev
Abstract: The problem of determining a spacewise-dependent right-hand side in a one-dimensional parabolic equation is considered. The method of conjugate gradients in combination with the method of finite differences with implicit time approximation with the weighting factor a € [0,1] is used for the numerical solution of the inverse initial-boundary value problem. The effectiveness of the proposed computational algorithm is confirmed by the results of the computational experiment for model problems with quasireal solutions, including problems with overdetermination conditions having random errors.
DOI: 10.25587/SVFU.2019.101.27249
Keywords: parabolic equation, inverse problem, source identification, difference scheme, conjugate gradient method.
REFERENCES
1. Samarskii A. A. and Vabishchevich P. N., Numerical Methods for Solving Inverse Problems of Mathematical Physics, De Gruyter (2007).
2. Kabanikhin S. I., Inverse and Ill-Posed Problems, Theory and Applications, De Gruyter (2011).
3. Lions J.-L. and Magenes E., Non-Homogeneous Boundary Value Problems and Applications, I, Springer, Berlin (1972).
4. Isakov V., Inverse Source Problems, Amer. Math. Soc., Providence, RI (1990) (Math. Surv. Monogr.; V. 34).
5. Vogel C. R., Computational Methods for Inverse Problems, SIAM (2002).
6. Cannon J. R., "Determination of an unknown heat source from overspecified boundary data," SIAM J. Numer. Anal., 5, No. 2, 275-286 (1968).
7. Isakov V., "Inverse parabolic problems with the final overdetermination," Commun. Pure Appl. Math., 44, No. 2, 185-209 (1991).
8. Solov'ev V. V., "Solvability of the inverse problems of finding a source, using overdetermination on the upper base for a parabolic equation," Differ. Equ., 25, 1114-1119 (1990).
9. Kostin A. B., "The inverse problem of recovering the source in a parabolic equation under a condition of nonlocal observation," Sb. Math., 204, 1391-1434 (2013).
10. Johansson B. T. and Lesnic D., "Determination of a spacewise dependent heat source," J. Comput. Appl. Math., 209, No. 1, 66-80 (2007).
11. Johansson B. T. and Lesnic D., "A variational method for identifying a spacewise-dependent heat source," IMA J. Appl. Math., 72, No. 6, 748-760 (2007).
The authors were supported by the grant from Russian Foundation for Basic Research (project No. 17-01-00732A).
© 2019 L. D. Su, V. I. Vasil'ev
12. Dhaeyer S., Johansson B. T., and Slodicka M., "Reconstruction of a spacewise-dependent heat source in a time-dependent heat diffusion process," IMA J. Appl. Math., 79, No. 1, 33—53 (2014).
13. Erdem A., Lesnic D., and Hasanov A., "Identification of a spacewise dependent heat source," Appl. Math. Model., 37, No. 24, 10231-10244 (2013).
14. Vabishchevich P. N., Additive Operator-Difference Schemes, Walter de Gruyter, Berlin; Boston (2014).
15. Liu J. and Tang J., "Variational iteration method for solving an inverse parabolic equation," Phys. Lett., A, 372, No. 20, 3569-3572 (2008).
16. Karchevsky A. L., "Reformulation of an inverse problem statement that reduces computational costs," Euras. J. Math. Comput. Appl., 1, No. 2, 4-20 (2013).
17. Samarskii A. A., Vabishchevich P. N., and Vasil'ev V. I., "Iterative solution of a retrospective inverse problem of heat conduction," Mat. Model., 9, No. 5, 119-127 (1997).
18. Vasil'ev V. I. and Kardashevsky A. M., "Iterative solution of the retrospective inverse problem for a parabolic equation using the conjugate gradient method," in: Rev. Sel. Pap. 6th Int. Conf. Numer. Anal. Its Appl. (Lozenets, Bulgaria, June 15-22, 2016), pp. 698-705, Springer, Cham (2017) (Lect. Notes Comput. Sci.; V. 1087).
19. Samarskii A. A., The Theory of Difference Schemes, Marcel Dekker, New York (2001).
20. Saad U., Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM (2003).
21. Vabishchevich P. N., "Iterative computational identification of a spacewise dependent the source in a parabolic equations," Inverse Probl. Sci. Eng., 1-23 (2016).
Submitted February 10, 2019 Revised February 26, 2019 Accepted March 1, 2019
Ling De Su and Vasilii I. Vasil'ev M. K. Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 42 Kulakovsky Street, Yakutsk 677000, Russia [email protected] and [email protected]
