где - смещенный прямоугольный импульс на шаг - кусочно-линейные функции на граничной окружности, причем и в остальных граничных узлах. В дальнейшем управления нумеруются единым индексом где - количество базисных управлений.
Восстановление плотности по рассчитанным решениям прямой задачи на множестве граничных узлов
проводилось по схеме, близкой к описанной выше схеме решения обратной задачи. При этом использовалась кусочно-постоянная модель плотности (в каждом треугольнике - константа).
Алгоритм решения обратной задачи в дискретной постановке состоял в следующем:
1.Рассчитывались сеточные "гармонические функции" где - количество граничных узлов:
где векторы - могут быть произвольными линейно-независимыми граничными векторами (их компоненты равны нулю во внутренних узлах),удовлетворяющих условию разрешимости задачи Неймана
2.Для каждой решалась задача граничного управления
(5.6)
относительно управлений . Как и в непрерывном случае, эти уравнения эквивалентны уравнениям, использующим только данные обратной задачи (состояния в обратной задаче неизвестны). Можно показать, что (5.6) выполняется тогда и только тогда, когда управление удовлетворяет линейному уравнению,
где билинейная форма явно выражается через данные обратной задачи
и условию
Во всех рассмотренных численных экспериментах задача граничного управления решалась с хорошей точностью.
3.Как и в случае непрерывной постановки BC - метод приводит к равенству
где правая часть явно вычисляется через данные обратной задачи. Подставляя сюда и заменяя на для нахождения значений плотности в каждом треугольнике получаем систему уравнений
Решая эту систему линейных алгебраических уравнений, мы восстанавливаем . Матрица системы оказывается плохо обусловленной, что естественно для многомерных обратных задач. Для ее решения использовался метод псевдообращения.
6. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ. РЕШЕНИЕ ДВУМЕРНОЙ ЗАДАЧИ
Для численного моделирования решения обратной задачи методом граничного управления был выбран математический пакет MatLab, так как в нем уже реализованы основные функции работы с матрицами, а также алгоритм триангуляции Делоне.
В качестве области была выбрана двумерная область – круг радиусом (рис. 6.1).
Для области построена оптимальная триангуляционная сетка со следующими параметрами: количество узлов сетки , количество граничных узлов , количество треугольников (рис. 6.2).
Построение оптимальной триангуляционной сетки Делоне обусловлено следующими параметрами: количеством узлов, количеством граничных узлов, количеством треугольников. Чем больше узлов сетке, тем лучше решается прямая задача, следовательно, улучшается решение задачи граничного управления. Чем больше граничных узлов, тем больше гармонических функций мы можем создать, следовательно, лучше будет решаться обратная задача – восстановление плотности. При увеличении количества треугольников сетки возрастает и количество неизвестных в системе уравнений (5.5), следовательно, результаты восстановления плотности станут хуже. Так как при «сгущении» сетки количество граничных узлов будет возрастать медленнее количества треугольников, то увеличение сетки приведет к ухудшению решения обратной задачи. В ходе численного моделирования была получена оптимальная сетка (рис.6.2).
Для решения поставленной задачи необходимо определить систему линейно-независимых управлений . В качестве были взяты функции следующего вида:
где - функция Рикера, смещенная на шаг (рис. 6.3), -количество шагов по времени, - характеристическая функция элемента границы (дуги s) на, изображенная на рис. 6.4. В проводимых численных экспериментах .
В ходе численного моделирования решения обратной задачи методом граничного управления были рассмотрены следующие примеры:
1.«Единичная плотность» ( рис. 6.5, 6.6).
Функция плотности: .
Задача нахождения массы. Масса области Относительная погрешность – 10-13.
Задача восстановления плотности. Относительная погрешность восстановления – 2,8%.
2.«Плотность с резким включением в центре» (рис. 6.7, 6.8).
Функция плотности: .
Масса области 0,0140. Относительная погрешность -10-12 .
Задача восстановления плотности. Относительная погрешность восстановления – 3,9%.
3. «Неодносвязная область с двумя большими резкими включениями» ( рис. 6.9, 6.10).
Функция плотности:
Масса области 0,0130. Относительная погрешность – 10-12.
Задача восстановления плотности. Относительная погрешность восстановления – 4,5%.
4. «Плотность, оцифрованная с монохромной картинки» (рис 6.11, 6.12).
Масса области 0,0163. Относительная погрешность -10-13 .
Задача восстановления плотности. Относительная погрешность восстановления – 3, 6%.
5. «Косинусоида» (рис. 6.13, 6.14).
Функция плотности: .
Масса области 0.0268. Относительная погрешность – 10-13.
Задача восстановления плотности. Относительная погрешность восстановления – 1,9%.
На основе написанного комплекса функций было создано учебное пособие «BC-method» для курса «Введение в теорию обратных задач» (Приложение 2). Основными задачами данного приложения являются:
1) Иллюстрация принципов работы метода граничного управления для решения обратной задачи параболического уравнения.
2) Возможность экспериментально показать зависимость работы метода от начальных условий (выбора сетки, плотности).
3) Возможность проведения различных численных экспериментов, используя удобный графический интерфейс для ввода входных параметров.
ЗАКЛЮЧЕНИЕ
В данной работе решалась двумерная обратная задача для параболического уравнения методом граничного управления. С помощью данного метода были проведены численные эксперименты для восстановления плотности и массы различных мембран. Погрешность при восстановлении коэффициентов в среднем 5%-6%, таким, образом, полученные результаты показали эффективность метода и возможность его использования при решении подобных задач.
Проведенные эксперименты позволяют предположить, что метод граничного управления применим и для ряда других обратных задач: обратные динамические задачи для системы уравнений теории упругости и системы Максвелла, имеющие приложения, в частности, в сейсморазведке и электроразведке, задач связанных с оптической томографией и других. Дипломный проект не является завершающим этапом проводимых исследований и планируется расширить круг решаемых задач.
На данный момент решаются только задачи небольшой размерностью (размерность зависит от количества управлений и количества узлов в триангуляционной сетке). Это ограничение связанно с недостатком вычислительных ресурсов. Поэтому, для того чтобы решать задачи больших размеров, необходимо использовать многопроцессорные системы. Для этого пишется программа на основе параллельных алгоритмов.
ЛИТЕРАТУРА
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. Segerlind L.J., Applied Finite Element Analysis. New York/ London/ Sydney/ Toronto: John Willey and Sons Inc, 1976. 392 p.
22. Стренг Г. Линейная алгебра и ее применения. М.: Мир, 1980. 454 с.
23. Бек ДЖ., Блакуэлл Б., Сент-Клэр Ч. .Некоторые обратные задачи теплопроводности.M.: Наука, 1989. 312 с.
24. Д. Норри, Ж. деФриз. Введение в метод конечных элементов.М.: Мир, 1981. 304 с.
25. Алифанов О.М. Обратные задачи теплообмена. М.: Машиностроение, 1988. 280 с.
26. Ануфриев И.Е., Смирнов А.Б., Смирнова Е. Н. MATLAB 7. СПб.: БХВ-Петербург, 2005. 1104 с.
ПРИЛОЖЕНИЕ 1
Алгоритм решения.
1. Вводим в консоли Matlab >>pdetool
2. В появившемся окне рисуем область и экспортируем данные о ней
(p, e, t).(рис П.1.1 и П.1.2).
Рис. П.1.1. Треангулированная область.
Рис. П.1.2. Экспорт параметров области.
3. [K,M]=assema(p,t,1,Ro,0); где Ro- плотность задается в виде функции или вектора.
4. Выполняем следующие команды:
K=full(K);
M=full(M);
[X,L]=eig(K,M);
U=matr_reakcii6(J,Q,p,e,X,L,40/256);
pinvP=pinv(U'*K*U);
pinvK=pinv(K);
Ct=U'*M*U;
Ck1=pinv(Ct)*Fi';
[Cpq,Qpq,C,E]=plotnosti(pinvK,pinvP,Ct,U,size(e,2),Ck1,p,e,t);
5. Вычисляем плотность
Rm=pinv(Qpq)*Cpq';
6.Прорисовка плотности.
pdeplot(p,e,t,'xydata',Rm','zdata',Rm','mesh','off');
7.Регуляризация
[Rm_1,nev]=reg(Qpq,Cpq,alfa,x0);%Rm_1 - плотность, nev - невязка
8. Просчет реальной плотности, если плотность задана ввиде функции
[email protected](x,y) eval(Ro);
for i=1:size(t,2)
Ro_real(i)=(R_fun(p(1,t(1,i)),p(2,t(1,i)))+R_fun(p(1,t(2,i)),p(2,t(2,i)))+R_fun(p(1,t(3,i)),p(2,t(3,i))))/3;
end;
9.Просчет погрешности востановления плотности
norma=norm(Rm-Ro_real',2)/norm(Ro_real,2);
Блок схема
Исходный код.
1. function U=matr_reakcii6(K,Q,p,e,X,L,T);
global k lj tt
L(1,1)=0;
L1=diag(L);
L2=L1.^(1/2);
shag_imp=T/Q;
G1=zeros(length(p),K);
for k2=0:1:K-1
k=k2;
[P,G,H,R]=assemb('boundmem64',p,e);
G1(:,k2+1)=full(G);
end
U=zeros(length(p),K*Q);
for impulse_at_time=0:1:Q-1
tt=shag_imp*(impulse_at_time+1);
L_resh=zeros(length(p),length(p));