Численное решение трехмерной обратной задачи для волнового уравнения методом граничного управления

Численное решение трехмерной обратной задачи для волнового уравнения методом граничного управления

...

Содержание
ВВЕДЕНИЕ 3
1. НЕФОРМАЛЬНАЯ ПОСТАНОВКА ЗАДАЧИ 8
2. ФОРМАЛЬНО-МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ 11
2.1. Прямая задача 11
2.2. Обратная задача 11
2.3. Метод граничного управления. Билинейные формы 12
2.4. Задача граничного управления 16
2.5. Задача нахождения массы 18
2.6. Задача восстановления плотности 19
3. ДИСКРЕТНАЯ ПОСТАНОВКА 21
3.1. Прямая задача. Метод конечных элементов (МКЭ) 21
3.2. Метод граничного управления в дискретном виде 23
3.3. Задача граничного управления в дискретном виде 23
3.4. Решение задачи нахождения массы в дискретном виде 25
3.5. Решение задачи восстановления плотности в дискретном виде 25
4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ РЕШЕНИЯ ТРЕХМЕРНОЙ ОБРАТНОЙ ЗАДАЧИ ДЛЯ ВОЛНОВОГО УРАВНЕНИЯ. 28
ЗАКЛЮЧЕНИЕ 36
ЛИТЕРАТУРА 38
ПРИЛОЖЕНИЕ 1. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ 40
П.1.1. Общее описание программы 40
П.1.2. Схема программы 43
П.1.3. Текст программы 49
ПРИЛОЖЕНИЕ 2. УЧЕБНОЕ ПРИЛОЖЕНИЕ «BC-METHOD» 64

ВВЕДЕНИЕ

Дифференциальные уравнения имеют большое прикладное значение, являясь мощным орудием исследования многих задач естествознания и техники: они широко используются в механике, астрономии, физике, во многих задачах химии, биологии. Это объясняется тем, что весьма часто законы, которым подчиняются те или иные процессы, записываются в форме дифференциальных уравнений, а сами эти уравнения, таким образом, являются средством для количественного выражения этих законов. Дифференциальные уравнения в частных производных представляют собой одну из наиболее сложных и одновременно интересных задач вычислительной математики. Эти уравнения характеризуются тем, что для их решения не существует единого универсального алгоритма, и большинство задач требует своего собственного особого подхода.
Все ди фференциальные уравнения с частными производными делятся на три основных типа: гиперболические, параболические и эллиптические.
Наиболее известным примером уравнений гиперболического типа является волновое уравнение. Данное уравнение моделирует малые поперечные колебания тонкой мембраны или струны, а также другие колебательные процессы в сплошных средах (акустика, преимущественно линейная: звук в газах, жидкостях и твёрдых телах) и электромагнетизме (электродинамике), а также находит применение и в других областях теоретической физики, например при описании гравитационных волн.
В различных областях науки и техники с целью познания закономерностей работы некоторого объекта или природного явления проводятся эксперименты самого различного вида. Цель эксперимента - выявление главных закономерностей явления и формирование на его основе некоторой математической модели. Очень часто на практике встречаются ситуации, когда объект исследования либо принципиально недоступен для наблюдения, либо проведение такого эксперимента дорого. Примерами таковых могут служить эксперименты по изучению внутреннего строения Земли, на основе которых можно было бы прогнозировать месторождения полезных ископаемых, предсказывать время и место разрушительных землетрясений. Отметим, что глубина самых глубоких шахт, пробуренных при помощи современнейшего оборудования, не превышает 13 км, а средний радиус Земли равен 6371 км. Таким образом, для непосредственных наблюдений колебаний Земли доступна лишь небольшая ее приповерхностная часть. При этом необходимо делать заключение о свойствах Земли (например, об изменении ее плотности с глубиной) по измеренным в ходе эксперимента косвенным проявлениям.
Другой пример - это медицинские исследования, направленные на выявление патологий внутренних органов человека. C открытием рентгеновских лучей человечество приобрело мощный инструмент исследования грудной клетки, костей, пищеварительного тракта, но в силу их неблагоприятного воздействия на ткани продолжались поиски менее вредного и более информативного способа изучения органов человека. Таким способом в настоящее время является ультразвуковое исследование (УЗИ), широко применяемое в медицинской практике и позволяющее достаточно просто выявлять патологии различных органов. В этом случае объект исследования также недоступен для непосредственного изучения. Мы судим о структуре и размерах органов лишь на основе косвенных данных измерений. В основе этого способа лежит анализ отраженных от органа волн.
У описанных выше примеров есть нечто общее - мы хотим определить причины, если известны полученные в результате экспериментов или наблюдений следствия. С точки зрения соотношения причина-следствие все задачи математического моделирования можно условно разделить на два больших класса: прямые задачи (известны причины, необходимо найти следствия) и обратные (известны следствия, нужно найти причины). К прямым задачам относятся, например, задачи расчета механических, тепловых, электромагнитных полей для тел, свойства которых и конфигурация известны. Эти задачи к настоящему времени достаточно хорошо изучены.
К обратным задачам относят задачи определения некоторых физических свойств объектов, таких, как плотность, коэффициент теплопроводности, упругие модули в зависимости от координат или в виде функций других параметров. Процедура решения таких задач, состоящих в обращении причинно-следственных связей, связана с преодолением серьезных математических трудностей. Успех ее сильно зависит как от качества и количества полученной из эксперимента информации, так и от способа ее обработки.

Группой студентов в состав, которой входит автор данной работы, было создано учебное пособие «BC-method» для курса «Введение в теорию обратных задач» (Приложение 2). Приложение «BC-method» решает следующие задачи:
1) Проведение численных экспериментов решения двумерной обратной задачи для волнового уравнения методом граничного управления.
2) Проведение численных экспериментов решения двумерной обратной задачи для уравнения теплопроводности методом граничного управления.
3) Иллюстрация рассчитанных примеров для волнового уравнения (двумерный, трехмерный случай) и для уравнения теплопроводности.
Численные эксперименты решения трехмерной обратной задачи для волнового уравнения, представлены в приложении «BC-method» в качестве примеров. Это связанно с тем, что для решения трехмерной задачинеобходимы мощные вычислительные ресурсы (оперативная память).
ЗАКЛЮЧЕНИЕ
В дипломном проекте решалась трехмерная обратная задача для волнового уравнения методом граничного управления, и проводились численные эксперименты по восстановлению плотности и массы. Полученные результаты показали эффективность метода и возможность его использования при обработке данных ультразвуковой томографии.
Проведенные эксперименты позволяют предположить, что метод граничного управления будет «работать» (то есть решать обратную задачу с достаточной точностью) и для ряда других обратных задач: обратные динамические задачи для системы уравнений теории упругости и системы Максвелла, имеющие приложения, в частности, в сейсморазведке и электроразведке. Дипломный проект не является завершающим этапом проводимых исследований и планируется расширить круг задач.
На данный момент решаются только «маленькие» задачи (расчетная область – куб со стороной , количество узлов триангуляционной сетки - 1331). Это ограничение связано с недостатком вычислительных ресурсов. В связи с этим, для того чтобы решать более крупные задачи, необходимо распараллелить производимые вычисления. В данный момент пишется программа для расчета «больших» задач на многопроцессорных системах.
По результатам проведенных исследований приняты к публикации статьи в журналах «Вестник НГУ» (г. Новосибирск), «Inverse Problems and Imaging» (American Institute of Mathematical Sciences) и опубликована в материалах IV научно-практической конференции «Обратные задачи и информационные технологии рационального природопользования» (г.Ханты-Мансийск). В 2008 году работа «Как услышать массу мембраны» была удостоена диплома I степени в подсекции "Математическое моделирование" XLVI Международной научной студенческой конференции "Студент и научно-технический прогресс" (г. Новосибирск). В 2009 г. работа «Численное восстановление плотности в обратной задаче для волнового уравнения методом граничного управления» заняла II место на XLVII Международной научной студенческой конференции "Студент и научно-технический прогресс".
ЛИТЕРАТУРА
1. Романов В. Г. Обратные задачи математической физики. М.: Наука, 1984. 263 с.
2. Белишев М. И., Благовещенский А. С. Динамические обратные задачи теории волн. СПб: Изд-во Санкт-Петербургского университета, 1999. 265 с.
3. Kabanikhin S. I. Definitions and examples of inverse and ill-posed problems // Journal of inverse and ill-posed problems. 2008. 16. №4. P. 317-357.
4. Гельфанд И. М., Левитан Б. М. Об определении дифференциального уравнения по его спектральной функции // Изв. АН СССР. Сер. мат. 1951. 15. № 4. С. 309-360.
5. Марченко В. А. Операторы Штурма-Лиувилля и их приложения. Киев: Наукова думка, 1987. 332 с.
6. Крейн М. Г. Об одном методе эффективного решения обратной краевой задачи // Докл. АН СССР. 1954. 94. № 6. С.767-770.
7. Благовещенский А. С. О локальном методе решения нестационарной обратной задачи для неоднородной струны // Тр. мат. ин-та АН СССР. 1971. 65. С. 28-38.
8. Романов В. Г. Устойчивость в обратных задачах. М.: Научный мир, 2004. 304 с.
9. Natterer F, Wubbeling F. A propagation-backpropagation method for ultrasound tomography // Inverse Problems. 1995. 11. P. 1225-1232.
10. Natterer F. Reectors in wave equation imaging // Wave Motion. 2008. 45. P. 776-784.
11. Кабанихин С. И. Проекционно-разностные методы определения коэффициентов гиперболических уравнений. Новосибирск: Наука (Сиб. отд.), 1988. 166 с.
12. Beilina L. Adaptive element/difference method for inverse lastic scattering waves // Appl. Comput. Math. 2002. 1. №2. P. 158-174.
13. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1974. 222 с.
14. Klibanov M. V., Timonov A. Carleman estimates for coefficient inverse problems and numerical applications. Utrecht (The Netherlands): VSP, 2004. 279 p.
15. Clason C., Klibanov M.V. The quasi-reversibility method for thermoacoustic tomography in a heterogeneous medium // SIAM J. SCI. COMPUT. 2007. 30. №2. P. 1-23.
16. Belishev M.I. Boundary control in reconstruction of manifolds and metrics (the BC-method) // Inverse Problems. 1997. 13. P. 1-45.
17. Belishev M.I. Recent progress in the boundary control method // Inverse Problems. 2007. 23. №5. P. 1-67.
18. Pestov L. N. On reconstruction of the speed of sound from a part of boundary // Journal of inverse and ill-posed problems. 1999. 7. №5. P. 481-486.
19. Lasiecka I., Lions J-L., Triggiani R. Non homogeneous boundary value problems for second order hyperbolic operators // J. Math. Pures Appl. 1986. 65. P. 149-192.
20. Bardos C., Lebeau G., Rauch J. Sharp sufficient conditions for the observation control and stabilization of the waves from the boundary // SIAM J. Contr. Opt. 1992. 30. P. 1024-1065.
21. kdjfdf
22. Segerlind L.J., Applied Finite Element Analysis. New York/ London/ Sydney/ Toronto: John Willey and Sons Inc, 1976. 392 p.
23. Стренг Г. Линейная алгебра и ее применения. М.: Мир, 1980. 454 с.
ПРИЛОЖЕНИЕ 1
ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
П.1.1. Общее описание программы
Для решения трехмерной обратной задачи для волнового уравнения методом граничного управления был написан комплекс функций, позволяющий проводить численные эксперименты решения задачи для различных входных параметров.
Для разработки данного комплекса была выбрана математическая среда MatLab. Все функции написаны на m-языке.
Структура комплекса функций (рис. П.1.1):
1) Main.m – главная функция, которая подготавливает данные для передачи дочерним функциям: Matrix_of_reaction.m, Ct_matrix.m, Density.m.
2) Area_KUB.m – функция, в которой задается сетка области (набор точек, создающий разбиение куба на более мелкие кубы).
3) Mat_e.m – функция, которая создает матрицу с номерами граничных треугольников.
4) Assema_3d.m – функция, которая рассчитывает матрицу масс и матрицу жесткости, а также плотность в элементах сетки (тетраэдрах).
5) Matrix_of_reaction.m – функция, которая решает прямую задачу и создает матрицу реакции - решение прямой задачи в граничных узлах для всех управлений в финальный момент времени.
6) Assemb_3d.m – функция, которая рассчитывает граничный вектор.
7) L_int1v2.m, L_intv2.m – функция, которая рассчитывает временную составляющую решения прямой задачи.
8) Ct_matrix.m – функция, которая рассчитывает матрицу скалярных произведений для базисных управлений.
9) Integral.m – функция, которая рассчитывает интеграл по времени, необходимый для расчета матрицы скалярного произведения.
10) Ricker.m – функция, которая считает функцию Рикера.
11) Int_t.m – функция, необходимая для расчета интеграла.
12) Density.m – функция, в которой производится расчет вектора-плотности.
13) Mass_ij.m – функция, рассчитывающая локальную матрицу масс.
14) Pdepir.m – функция, рассчитывающая объем элемента сетки (тетраэдра).
15) Drawing – функция, отрисовывающая плотность области на срезе плоскостью z=0.5.
Рис. П.1.1. Структура комплекса функций.
Входные параметры:
- радиус расчетной области (вещественное число).
- плотность области, заданная в виде аналитической функции (строковая переменная).
- количество управлений по пространству (натуральное число).
- количество управлений по времени (натуральное число).
- финальный момент времени (вещественное число).
Выходные параметры:
- масса области (вещественное число).
- относительная погрешность расчета массы области (вещественное число).
- восстановленная плотность (вектор размерностью ).
- относительная погрешность восстановления плотности области (вещественное число).
Для вызова программы необходимо:
1) Поместить исполняемые файлы в рабочий каталог (по умолчанию «..\MATLAB\R2006b\work»).
2) Запустить среду Matlab (выполняемый файл matlab.exe).
3) Из командной строки среды MatLab вызвать функцию Main с заданными входными параметрами (>>[mas,delta_mas,Ro_restore,delta_Ro_restore]=Main(Ro,S,Q,T)).
Размерность решаемых задач зависит от количества узлов сетки и количества управлений. Чем больше размерность поставленной задачи, тем больше вычислительных ресурсов требуется для ее решения. В связи с тем, что трехмерная задача имеет большую размерность (работа с матрицами порядка 7000х7000), вычислительных ресурсов ПК не хватает для расчета данной задачи. Таким образом, для расчетов использовался суперкомпьютер SunFire 15K.
П.1.2. Схема программы
П.1.3. Текст программы
1) Main.m
function [mas,delta_mas,Ro_restore,delta_Ro_restore]=Main(Ro,S,Q,T)
%задаем набор точек, создающий разбиение куба на более мелкие кубы
[p]=area_KUB;
%для данного набора точек создаем триангуляционную сетку
t=delaunay3(p(1,:), p(2,:), p(3,:));
t=t';
%записываем номера граничных треугольников
e=mat_e(p,t);
%расчитываем матрицы масс и жесткости, а также находим реальное значение
%плотности в элементах триангуляционной сетки (тетраэдрах)
[K,M,Ro_real]=assema_3d(p,t,Ro);
K=full(K);
M=full(M);
%решаем задачу на собственные значения
[X,L]=eig(K,M);
%решаем прямую задачу - находим матрицу реакции
U=Matrix_of_reaction(S,Q,p,e,X,L,T);
%рассчитываем матрицу потенциальной билинейной формы
P=U'*K*U;
%рассчитываем матрицу скалярного произведения волн
[Ct]=Matrix_of_Ct(S,Q,p,e,X,L,T);
%находим управление "нагоняющее" единицу
P=(P+P')/2;
[x,l]=eig(P);
[a,b]=min(abs(diag(l)));
Ck1=x(:,b)/mean(U*x(:,b));
%восстанавливаем массу области
mas=Ck1'*Ct*Ck1;
%находим относительную ошибку восстановления массы
delta_mas=norm(mas-sum(sum(M)),2)/norm(sum(sum(M)),2);
%находим псевдообращение матрицы жесткости, матрицы потенциальной
%билинейной формы
pinvK=pinv(K);
pinvP=pinv(P);
%находим матрицы необходимые для нахождения плотности
[Qpq,Cpq,C,E]=Density(pinvK,pinvP,Ct,U,120,Ck1,p,e,t);
%находим псевдообращение матрицы Qpq
pinvQpq=pinv(Qpq);
%восстанавливаем плотность
Ro_restore=pinvQpq*Cpq';
%находим относительную ошибку восстановления плотности
delta_Ro_restore=norm(Ro_restore-Ro_real,2)/norm(Ro_real,2);
%отрисовываем результат восстановления
Drawing(p,t,Ro_restore);
2) Area_KUB.m
%задает сетку области
function [p]=area_KUB
kol=0;
for i=0:0.1:1
for j=0:0.1:1
for k=0:0.1:1
kol=kol+1; p(1,kol)=i; p(2,kol)=j; p(3,kol)=k;
end
end
end
3) Mat_e.m
% определяет номера граничных треугольников
function y=mat_e(p,t);
e=[]; kol=0; m=max(max(p));
for tet=1:1:size(t,2)
k=0;
for i=1:1:4
if p(1,t(i,tet))==0
k=k+1;
end
end
if k==3
kol=kol+1;
k=0;
for i=1:1:4
if p(1,t(i,tet))==0
k=k+1;
e(k,kol)=t(i,tet);
end
end
end
end
for tet=1:1:size(t,2)
k=0;
for i=1:1:4
if p(1,t(i,tet))==m
k=k+1;
end
end
if k==3
kol=kol+1;
k=0;
for i=1:1:4
if p(1,t(i,tet))==m
k=k+1;
e(k,kol)=t(i,tet);
end
end
end
end
for tet=1:1:size(t,2)
k=0;
for i=1:1:4
if p(2,t(i,tet))==0
k=k+1;
end
end
if k==3
kol=kol+1;
k=0;
for i=1:1:4
if p(2,t(i,tet))==0
k=k+1;
e(k,kol)=t(i,tet);
end
end
end
end
for tet=1:1:size(t,2)
k=0;
for i=1:1:4
if p(2,t(i,tet))==m
k=k+1;
end
end
if k==3
kol=kol+1;
k=0;
for i=1:1:4
if p(2,t(i,tet))==m
k=k+1;
e(k,kol)=t(i,tet);
end
end
end
end
for tet=1:1:size(t,2)
k=0;
for i=1:1:4
if p(3,t(i,tet))==0
k=k+1;
end
end
if k==3
kol=kol+1;
k=0;
for i=1:1:4
if p(3,t(i,tet))==0
k=k+1;
e(k,kol)=t(i,tet);
end
end
end
end
for tet=1:1:size(t,2)
k=0;
for i=1:1:4
if p(3,t(i,tet))==m
k=k+1;
end
end
if k==3
kol=kol+1;
k=0;
for i=1:1:4
if p(3,t(i,tet))==m
k=k+1;
e(k,kol)=t(i,tet);
end
end
end
end
y=e;
4) Assema_3d.m
%расчет матрицы масс, матрицы жесткости, вектора плотности (которую необходимо восстановить) в элементах сетки.
function [K,M,Ro_real]=assema_3d(p,t,Ro);
K=sparse(size(p,2),size(p,2));
M=sparse(size(p,2),size(p,2));
[email protected](x,y,z) eval(Ro);
for tet=1:1:size(t,2)
A=[1 p(1,t(1,tet)) p(2,t(1,tet)) p(3,t(1,tet));...
1 p(1,t(2,tet)) p(2,t(2,tet)) p(3,t(2,tet));...
1 p(1,t(3,tet)) p(2,t(3,tet)) p(3,t(3,tet));...
1 p(1,t(4,tet)) p(2,t(4,tet)) p(3,t(4,tet))];
vol=abs(det(A))/6;
Amin=A^(-1);
C=sum(A(:,2:4))/4;