Содержание к диссертации
Введение
Глава 1. Оценка производных по параметрам математических ожиданий функционалов от диффузионных процессов с условием поглощения на границе 16
1.1. Постановка задачи 16
1.2. Применение метода Эйлера для получения оценок функционалов от диффузионных процессов 19
1.3. Оценка производных ди/дв 22
1.4. Метод FBVP1 для определения оценок производных BHJ\& ди/дв 24
1.5. Моделирование траекторий процессов X и Л 51
1.6. Численный пример 53
1.7. Метод FBVP2 для определения оценок производных вир а ди/дв 56
1.8. Анализ погрешности аппроксимации производных по параметрам 63
1.9. Численная проверка метода FBVP2 74
Глава 2. Оценка производных по параметрам математических ожиданий функционалов от диффузионных процессов с условием отражения на границе 78
2.1. Моделирование диффузионного процесса в области с отражающей границей 81
2.2. Дифференцирование по параметру решения уравнения Скорохода 85
2.3. Оценки производных ди/дв 96
2.4. Решение модельной задачи 100
Глава 3. Решение обратных задач для уравнения теплопроводности с использованием оценок ди/дв 103
3.1. Определение неизвестных коэффициентов в краевой задаче для уравнения теплопроводности 103
3.2. Решение модельной задачи 111
Глава 4. Минимизации дисперсии функционала от диффузионного процесса с использованием преобразования краевой задачи для параболического уравнения 115
4.1. Вариант метода для поглощающей границы 115
4.1.1. Оценка математического ожидания функционала и дисперсия оценки
4.1.2. Преобразование краевой задачи 124
4.1.3. Численные эксперименты 127
4.2. Вариант метода для условия отражения на границе 134
4.2.1. Постановка задачи минимизации дисперсии оценки 134
4.2.2. Преобразование краевой задачи 140
4.2.3. Численные эксперименты 143
Глава 5. Оценка решения линейного уравнения теплопроводности с разрыв ными коэффициентами 147
5.1. Аппроксимация обобщенного решения уравнения теплопроводности
с разрывными коэффициентами 147
5.2. Сходимость приближенного решения краевой задачи и его статистическая оценка 152
5.3. Расчет теплового состояния сотовой теплозащитной панели
5.3.1 Расчет модельного примера 161
5.3.2 Расчет температуры сотовой панели по данным физического эксперимента 163 Заключение 167
Список работ, в которых отражены результаты диссертации 168
Список литературы
- Метод FBVP2 для определения оценок производных вир а ди/дв
- Дифференцирование по параметру решения уравнения Скорохода
- Решение модельной задачи
- Вариант метода для условия отражения на границе
Введение к работе
Актуальность темы исследования. Задачи эволюционного типа со случайными параметрами, математическое описание которых содержит случайные процессы диффузионного типа, встречаются во многих областях науки и техники. К таким задачам относятся, например, задачи автоматического управления, теории надежности, финансовой математики и др. Их решение обычно сводится к оценке требуемых усредненных характеристик, например, вероятностей каких-либо событий, оценок моментов распределений функционалов от случайных процессов. Как правило, успешное решение такого типа задач возможно с использованием метода статистического моделирования на основе численного решения стохастических дифференциальных уравнений (СДУ).
Известны также применения метода Монте-Карло с использованием решения СДУ и для нахождения решений детерминированных задач таких как, например, линейных эллиптических и параболических уравнений на основе вероятностных представлений их решений. Эти вероятностные представления выражаются в виде математических ожиданий функционалов от решений СДУ (Дынкин Е.Б. Марковские процессы. М.: Физматгиз, 1963. 860с). Применение метода Монте-Карло в к решению детерминированных задач целесообразно при их большой размерности, т.к. для нахождения решения не требуется построение сетки по пространственным переменным, и с ростом размерности вычислительные затраты растут незначительно. Параболические уравнения большой размерности возникают, например, в задачах финансовой математики (Ширяев А.И. Основы стохастической финансовой математики. Т.1. М.: Фазис, 1998. 489с.) и в задачах диффузионной динамики наночастиц и макромолекул (Mikkelsen A. et al. Brownian dynamics symulation of needle spring chains. Physica A. 1998. 253. P. 66-76.).
Применение метода Монте-Карло может быть также полезным, когда требуется найти решение и его производные по параметрам или по пространственным переменным в одной или нескольких точках области, а для решения сеточным методом, в силу особенности решаемой задачи, требуется очень большое количество узлов. Существенным аргументом в пользу методов Монте-Карло является их относительная простота реализации и
хорошая эффективность распараллеливания.
Важнейшей характеристикой практически любой математической модели является ее параметрическая чувствительность. Для ее исследования возможны как экспериментальные методы на основе физических экспериментов, так и математические, требующие определения производных по параметрам решений уравнений модели. При этом, как правило, математические методы менее дорогостоящие по сравнению с экспериментальными, но они требуют дополнительных теоретических исследований и построения алгоритмов расчета.
В задачах, в которых основную роль играют случайные процессы, определение параметрической чувствительности сводится к нахождению математических ожиданий производных по параметрам функционалов от этих случайных процессов. При этом, если по условиям задачи требуется, чтобы случайный процесс находился в заданной в области, важное значение имеют заданные граничные условия. Так, например, если на границе задано условие поглощения, то при вычислении производных следует учитывать производную по параметрам времени первого выхода диффузионного процесса из области. Если на границе задано условие отражения, то следует учитывать зависимость от параметров локального времени пребывания диффузионного процесса на границе.
Основные цели диссертационной работы:
Исследование возможности дифференцирования по параметрам математических ожиданий функционалов от диффузионных процессов в областях с поглощающими и отражающими границами.
Разработка методов на основе численного решения СДУ оценки математических ожиданий функционалов от диффузионных процессов и их производных по параметрам.
Методы исследования. В данной работе использовались методы теории вероятностей, методы численного моделирования и анализа стохастических дифференциальных уравнений, методы оптимизации, методы функционального анализа, методы теории решения обратных задач, методы исследования краевых задач для параболических уравнений.
Научная новизна. Все основные результаты диссертации являются но-
выми. В работе впервые найден способ вычисления производных по параметрам математических ожиданий функционалов от диффузионных процессов, содержащих время первого выхода процесса из области, на основе которого построены два метода определения оценок производных по параметрам математических ожиданий функционалов от диффузионных процессов в областях с поглощающими границами. Исследована возможность дифференцирования по параметрам функционалов от диффузионных процессов с условием отражения на границе и разработан численный метод для получения оценок производных по параметрам математических ожиданий таких функционалов. Получена формула для предельного значения дисперсии оценки математического ожидания функционала от диффузионного процесса при убывании длины шага в методе Эйлера. Разработан метод уменьшения дисперсии функционалов от диффузионных процессов, использующий преобразование краевой задачи для параболического уравнения. Разработан метод для получения оценок решения параболического уравнения с разрывными коэффициентами, основанный на численном решении СДУ.
Теоретическая ценность и практическая значимость. В диссертации показана возможность дифференцирования по параметрам математических ожиданий функционалов от диффузионных процессов в областях с условиями поглощения и отражения на границе. В случае поглощающих границ получены формулы для для определения производных по параметрам математических ожиданий функционалов, содержащих время первого выхода диффузионного процесса из области. Следует отметить, что эти формулы не содержат производные по параметрам от времени первого выхода процесса из области. Основной проблемой, которая возникает при дифференцировании математических ожиданий функционалов в случае отражающих границ, является необходимость дифференцирования по параметрам локального времени пребывания диффузионного процесса на границе. Автором показано, что при определенных условиях на входные данные задачи такое дифференцирование возможно. На этой основе был разработан метод для определения производных функционалов от диффузионных процессов в случае, когда на границе задано условие отражения.
В диссертации на основе полученной формулы для дисперсии функционалов от диффузионных процессов построен метод минимизации дисперсии таких функционалов.
Разработанные автором методы определения производных по параметрам могут применяться для решения задач параметрической оптимизации градиентными методами. В работе такая оптимизация используется при решении обратной задачи для уравнения теплопроводности и в задаче минимизации дисперсии.
В диссертации разработан метод на основе численного решения СДУ, с помощью которого можно определять решения краевых задач для параболических уравнений с разрывными коэффициентами в заданных точках области. Данный метод используется на практике для определения теплового состояния теплозащитных покрытий фюзеляжей самолетов, содержащих теплозащитные панели сотового типа.
Основные результаты, выносимые на защиту:
-
Получены две формулы для вычисления производных по параметрам математических ожиданий функционалов от диффузионных процессов, содержащих время первого выхода случайного процесса из области. На основе этих формул построены два метода для получения оценок производных по параметрам математических ожиданий функционалов от диффузионных процессов с условием поглощения на границе.
-
Построен метод оценки производных по параметрам математических ожиданий функционалов от диффузионных процессов в областях с условием отражения на границе.
-
Получена формула для предельного значения дисперсии оценки математического ожидания функционала от диффузионного процесса при убывании длины шага в методе Эйлера.
-
Разработан метод минимизации дисперсии оценки математического ожидания функционала от диффузионного процесса, основанный на преобразовании краевой задачи для параболического уравнения.
-
Разработан метод оценки решений краевых задач для уравнения теплопроводности с разрывными коэффициентами.
-
Решена задача оценки теплового состояния сотовых конструкций фю-
зеляжа самолета на основе численного решения СДУ.
Апробация работы. Основные результаты диссертационной работы неоднократно докладывались и обсуждались на семинаре "Методы Монте-Карло в вычислительной математике и математической физике "под руководством чл.-корр. РАН ГА. Михайлова в ИВМиМГ СОРАН, докладывались на семинаре кафедры теории вероятностей МАИ под руководством д.ф.-м.н., проф. А.И. Кибзуна, докладывались на четвертом , пятом и шестом Санкт-Петербургских семинарах по стохастическому моделированию (Санкт Петербург, 2001, 2005, 2009 гг.); на международной конференции по вычислительным методам "International Conference on Computational Methods in Science and Engeneering CMMSE -(CMMSE-2002)"(Аликанте, Испания, 2002 г.); на четвертом международном семинаре по методам Монте-Карло "IVth IMACS Seminar on Monte Carlo methods September 15-19, 2003, Berlin"(Берлин, Германия, 2003 г.); на международной конференции, посвященной 100-летию академика И.И. Векуа "Дифференциальные уравнения, теория функций и приложения"(Новосибирск, 2007 г.); на международной конференции, посвященной 100-летию С.Л. Соболева "Дифференциальные уравнения, функциональные пространства, теория приближений"(Новосибирск, 2008 г.); на восьмом Всемирном конгрессе по вычислительной механике "8th World Congress on Computational Mechanics WCCM8"h пятом Европейском конгрессе по вычислительным методам в прикладных науках и технике "5th European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2008"(Венеция, Италия, 2008 г.); на пятой международной конференции по машиностроению и аэрокосмической технике "The 5th International Conference on Mechanical and Aerospace Engineering (ICMAE 2014)"(Мадрид, Испания, 2014 г.); на международных конференциях "Актуальные проблемы вычислительной и прикладной математики": АПВМ2014, АПВМ2015 (Новосибирск, 2014, 2015 гг.); на 19-й конференции по системотехнике и компьютерным технологиям "19th International Conference on Circuits, Systems, Communications and Computers, (CSCC 2015)" (Закинтос, Греция, 2015 г.); по прикладным методам статистического анализа "AMSA 2015"(Новосибирск, 2015 г.).
Публикации. По теме диссертации опубликовано 21 работ. Из них 12 работ опубликовано в журналах, рекомендованных ВАК. Список работ автора по теме диссертации приведен в конце автореферата.
Структура и объем диссертации. Диссертация состоит из введения, 5 глав, заключения, списка литературы и приложения. Общий объем рпботы составляет 188 страниц.
Метод FBVP2 для определения оценок производных вир а ди/дв
В работах [79] - [82] показано, что погрешность аппроксимации метода Эйлера при вычислении математических ожиданий диффузионных процессов с поглощением на границе есть величина порядка 0{yh). Как было установлено в многочисленных исследованиях (см., например, [68] - [72], [79]), основная причина такой низкой точности состоит в следующем. Для любых двух последовательных значений Х{, ХІ+І, полученных на одном шаге метода Эйлера, принадлежащих G, существует ненулевая вероятность того, что соответствующий непрерывный процесс Xt на отрезке [ , +i], определяемый уравнением Xt = X% + ha(ti,Xi)(t - U) + Y AU,Xi){W(t) - W(U))j (1.10) 3=1 со значениями на границах отрезка Х{: Х{+\, есть процесс типа броуновского моста, который на интервале ( , +i) может выйти из области G и снова в нее вернуться. То есть численный метод может продолжать работу внутри области после выхода реального процесса из области.
Если известна вероятность выхода случайного процесса (1.10) из области на г-м шаге р(і, і +1) или ее достаточно хорошая оценка, то в целях улучшения сходимости в работах [79] - [82], [100] предлагается моделировать выход процесса (1.10) из области. В [70] -[72] (Buchmann F.M., Petersen W.P.) путем вычислений на простых примерах для одномерного случая показано, что порядок метода Эйлера может быть равен единице, если на г-м шаге (і = 0,... ,TN Л N - 1 моделировать случайную величину , распределенную по закону Бернулли с вероятностью "успеха" p(i,i + 1). При этом на шаге і возможны следующие ситуации: а) Х{ Є G, Xi+1 G, в этом случае ясно, что выход из области произошел; б) ХІ,ХІ+1 Є G, и, если при моделировании оказывается = 0, то в этом случае осуществляется переход к следующему шагу; если = 1, то фиксируется выход случайного процесса из области и траектория заканчивается. В работе [100] для случая, когда область G является полупространством, соответственно при этом dG - полуплоскость, получена формула для вероятности р(г,і + 1) = ехр ( -2 -- " ) (І-"] (п,X, - z)(n,Xi+1 - z) nT(7(7T(Xi)n/l где п - вектор нормали dG; z - проекция Х{ на dG.
В работе [80] для случая гладкой границы, гладкости не ниже С5, формула (1.11) используется в качестве оценки вероятности р(і, і + 1) в некоторой заданной окрестности границы, в которой область аппроксимируется полупространством, а граница - плоскостью. Там же утверждается, что для точек вне этой окрестности выход случайного процесса из области маловероятен, и дана оценка этой вероятности. Для варианта метода, в котором производится аппроксимации границы плоскостью с целью моделирования выхода процесса из области с использованием вероятности р(і, і + 1) в [80] доказан первый порядок сходимости.
В дальнейшем, для сокращения записи формул, будем считать, что в уравнении (1.1) присутствует один скалярный параметр в, отметив при этом, что все результаты легко переносятся на случай нескольких параметров. Также иногда аргумент в в списках аргументов функций будет опускаться. Если формально продифференцировать по в равенство (1.2), то получим + ( (т,хт,в)гт + (т,хт,в) )Хт т ТЛт + / ( (v,Xv,0)Zv + (v,Xv,0)yv\ + Ф(0) , (1.12) t где т(в + Щ-т(в) ,.., v п. дф Ф(в) :=tlEt + - (f(r 0) + f (т в))Хт т) (1.13) если этот предел существует), и где процесс Z.W = / (! ( ) + ) + / ( ( ) + %) mv) (1.14) является среднеквадратической (с.к.) производной - по параметру решения СДУ (1.1). Известно (см., например, [17]), что предположения А -Ев разделе 1.1 обеспечивают существование с.к. производной Zs{6) = -g -, которая может быть получена в результате решения системы СДУ (1.1), (1-14).
Для определения первого слагаемого в (1.12) требуется только знать процесс (Xm(9)}Zm(9))} который является марковским, т.к. система СДУ (1.1), (1-14) удовлетворяет условиям существования сильного решения.
При определении второго слагаемого в (1.12) возникает необходимость учитывать зависимость г от в. Доказательство существования предела в (1.13) и методика определения Ф(в) даны в последующих разделах. 1.4 Метод FBVP1 для определения оценок производных вида ди/дв
Ниже дано описание метода FBVP11 (см. [86], [87]) для определения статистических оценок производных математических ожиданий функционалов от диффузионных процессов в областях с поглощающими границами, основанного на применении формулы Ито к некоторой достаточно гладкой функции д(х). При этом функция д должна быть равной нулю на границе области 9G, а функции Lg ( L - дифференциальный оператор (1.6)) и dtGij не равны нулю одновре-менно при любых (t, ж, 9) Є ST Х U.
Оценки математических ожиданий функционалов, содержащих производные дт/89 получаются с помощью уравнения, полученного в результате применения формулы Ито к функции д. Поскольку формула Ито содержит недифференци-русмый интеграл Ито по винеровскому процессу, в верхнем пределе которого присутствует т, то предлагается аппроксимировать интеграл Ито с.к. интегралом, в котором белый шум аппроксимируется стационарным гауссовским процессом с экспоненциальной корреляционной функцией, согласованной с шагом интегрирования в методе Эйлера. Идея аппроксимации белого шума гауссовским процессом рассматривалась в [60]. При этом брался стационарный гауссовский процесс Л с корреляционной функцией R\{ti,t2) = D exp(-/3\t2 - h\). (1.15) Согласно [1], [60] стационарный гауссовский процесс Л с нулевым математическим ожиданием и экспоненциальной корреляционной функцией может быть 1FBVP - First Boundary Value Problem
Дифференцирование по параметру решения уравнения Скорохода
Для численного моделирования траекторий системы СДУ (2.1), (2.2) предлагается использовать вариант схемы Эйлера с ортогональным проектированием численного решения на G. В работах [69], [74] была исследована погрешность оценки решения параболического уравнения без источника, получаемая на основе этого метода, и был установлен порядок аппроксимации 0{h/ ). При этом рассматривался другой вид функционала, т.е. такой, в котором отсутствует интеграл от источника.
В работе [69] был также исследован так называемый симметризованный метод, в котором при выходе траектории из области на каком-либо шаге численного моделирования, очередная точка на данном шаге берется внутри области симметрично относительно границы. Для симметризованного метода установлен первый порядок сходимости. Но при этом остается открытым вопрос об оценке математических ожиданий функционалов, содержащих интегралы по локальному времени.
Пусть h - шаг интегрирования и положим N = (Т — t)/h. На отрезке [t,T] зададим сетку с узлами t{ = t + г/г, (г = 0,..., N). Свяжем с узлами этой сетки следующие обозначения: Х{, к{ — аппроксимация траектории системы СДУ (2.1)-(2.2) в точке ti\ Aj+1X, А{+1К — изменение приближенных значений X, Кш, за один шаг численного интегрирования; ді — значение g(ti,X{) для некоторой функции g(t,x); щ — единичный вектор внутренней нормали в точке p(X+1). Далее, для указания k-Q. компоненты вектора v будет использоваться обозначение (v)k .
Численное моделирование траекторий системы СДУ (2.1), (2.2) в узлах заданной сетки осуществляется по следующей схеме Х+1 =Хг + hat + Vh T&, (2.9) Xi+ X+ 1 + iAi+ m, (2.10) кг+1 = кг + Аг+1К, (2.11) где j - нормальные N(0,1) случайные величины; А{+1К = [ (JQ+ 1)]-. Заменив в функциях (2.4), (2.5) интегрирование суммированием по узлам сетки, определим соответствующие сеточные функции У , Zj, заданные в узлах (ti,Xi)(i = 0,...,N), 1, г = 0, У і = г-1 _ (2.12) exp(J2(hck + VkXdG(Xk)Ak+iK)), і = 1,...,7V, /г=0 О, г = 0, Z, = г-1 _ _ (2.13) Y.(hfk + lkXdG{Xk)Ak+iK)Yk) i = l,...,N. к=0 В качестве приближенного решения краевой задачи (2.6) — (2.8) принимается оценка uh(t,x) = Et (XN)YN + ZN), (2.14) которая получается в результате моделирования траекторий системы (2.1), (2.2) по формулам (2.9) — (2.13). В следующей теореме доказывается, что оценка (2.14) имеет порядок точности 0{Ш) . Теорема 6. Пусть выполнены условия А2; Б27 В2. Тогда при достаточно малых h существует константа С, зависящая от Т — t, а, о , f такая, что \uh(t,x)-u(t,x)\ Cht. (2.15) Доказательство. Рассмотрим разность uh{t, х) - u{t, х) = EtjX( p(XN)YN + ZN) - u{t, x) = EtjX(u{tN, XN)YN + ZN) - u{t, x) N-l = J2 { АФг+ъХг+1)г+1 + гг+1-и(и}х г-г ) . (2.16) Найдем оценку выражения Ei)X ( ui+1Yi+1 + Zl+1 - щУі - Zi ) . (2.17) Используя разложение щ+1 в точке (ti,X{) в ряд Тейлора до членов, содержащих производные второго порядка функции и: включительно и определение Y{, Zi7 получим /г Э2и; , А _ч N- d2Ui ICi +г]гХдс(хг)Аг+1К) + Z,+ (hfi + -/гХдс(хг)Аг+1К)Уг. (2.18) По построению, АІ+1К есть величина порядка 0{h/ ). Запишем равенство, которое получается в результате разложения экспоненты при достаточно малых h ехр(/гсг + г)гХдс{Хг)Аг+1К) = 1 + Нсг + г)гХдс{Хг)Аг+1К +\{Ъъ + т]гХдс(Х Аг+1К) 2 + O(Zii). (2.19) При получении оценки (2.15) мы используем тот факт, что при некотором ho О выполняется [112] sup EttX((kN)q) ос \/q Є N, (2.20) /i /io и, что вероятность выхода моделируемой траектории из области есть величина порядка 0{Ьл) [79] Р{Хг Є G,X+1 iG} Ch). (2.21) Из (2.18), (2.19) с учетом (2.20), (2.21) и того, что в точке (ti,Xi) выполняется уравнение (2.6), получим (f -E A.+1X dxk Е ,ж [Ui+1Уі+1 + Zi+1 — ЩУІ — Zi ) = Et, 1 d2Ui + 2 2 о о (A,+1X)fc(A,+1X)/jy, + YlClUlh + гг]гигХдс(Хі)Аг+1К к,І Хк Xl (Til I +Yifih + УгЪХдс(Х Аг+1К + YtJ2 - (Ьі+1Х гнхдсІХі)Аг+1К] +0(/i5) = Et, +ЪХдс(Х, диг дхк УЛг+1К Пі)кХк \й(х-1 ) + ViUiXdG (X i 2.22 + 0(h ). Ясно, что при Х{ Є G и Х+ 1 Є G, выражение под знаком математического ожидания в правой части (2.22) равно нулю, а при Х{ Є dG и Х+ 1 G оно обращается в ноль в силу выполнения граничного условия (2.8). В случае, когда Х{ Є G и Х+ 1 G это выражение равно 9и,: к к Но здесь следует отметить, что справедлива оценка N-1 диг - J2 (УД+1К р Іш +1 )) CE t k N )ti г=0 к Хк которая следует из (2.20) и (2.21). Из (2.22) и (2.23) следует, что \uh(t,x) -u(t,x)\ Chk. Теорема доказана. 2.2 Дифференцирование по параметру решения уравнения Скорохода 2.23) Рассмотрим сначала детерминированный случай. Пусть дана неслучайная функция r(t, в) с областью определения [0, Т] х [61, $2], принимающая значения в Rd, и г(0, в) Є О при всех значениях параметра в Є [#і, #2] Будем предполагать, что для функции г выполнено условие Д2) функция г непрерывна по t на отрезке [О, Т] и имеет производную по в в каждой точке интервала (#i, #2), непрерывную по і и по #. Решение уравнения Скорохода состоит в нахождении непрерывной функции /І : [О, Т] х [#і, 6Q\ — G И функции ограниченной вариации к : [О, Т] х [#1? #2] — Rd таких, что / (0, #) = 0 (і = 1,..., d) при всех в Є [ , ],и1г изменяется только когда /І Є dG. При этом для всех (t,6) Є [0,Т] х [#1,6 2] должны выполняться равенства: (i = r + k, (2.24) t fc( ,0) = f XdG(Ks,0))d\k(s,e)\, (2.25) о t h(t,0)= I n(n(s,0))d\k(s,0)\, (2.26) 0 где &(,#) при фиксированном в есть полная вариация функции k(t,6) на отрезке [0,]. По определению полной вариации ft(,#) = sirpj ft(j,#) — k(ti-i,0)\, к где супремум берется по всем разбиениям вида 0 = to t\ tn = t. Теоремы существования и единственности решений уравнений вида (2.24) даны в работах [102], [ПО], [116]. В работах [ПО], [116] доказательство существования и единственности решения уравнения (2.24) проводилось сначала для ступенчатых функции и затем осуществлялся переход к пределу. Для обоснования дифференцируемости по параметру решения уравнения (2.24) тоже будем сначала рассматривать ступенчатые функции. Зададим равномерное разбиение 0 = to t\ ім = Т отрезка [0,Т] с шагом h = ti — ti-i = l/2m, (і = 1,... ,7V — 1), где число узлов разбиения N = [2mT], если 2mT - [2mT] = 0 и N = [2тТ] + 1 в противном случае. Непрерывная функция r{t,9) на отрезке [0,Т] равномерно по в Є [#і,#2] может быть аппроксимирована ступенчатыми функциями, которые определяются по правилу: rm(t,Є) = г(к2 т}Є) при к2 т t {k + l)2"m, k = 0,...,N-l. Рассмотрим уравнение Скорохода для гто: Для функции гт решение уравнения (2.27) записывается явно [ПО] .MJrM t (2,8) 0(u (1=1 О) л. Аг ( — ) ) — І г-±± 0, 0 2І, М , 0)= кт{ ,Є) + р{ііт(#, в) + Агт( , в) ) (2.29) _і, (і=1 и) - Аг (— в) — t ГДЄ АЛ ГутДт т і ") m\QmJ т\2т "/ В дальнейшем, производную по параметру в некоторой функции а (в) будем обозначать как два. Из дифференцируемости по в функции г, дифференцируемое по х функции р(х) и формул (2.28), (2.29) ясно, что функции гто, кт дифференцируемы по в. В силу дифференцируемости р(х) ( проекции точки х на G) и функции г по в можно последовательно для каждого интервала разбиения отрезка [0,Т] продифференцировать по в равенства (2.28), (2.29). В результате получим производные функций гт, кт по параметру
Решение модельной задачи
Если вычислительные затраты, связанные с преобразованием краевой задачи, существенно меньше вычислительных затрат, требуемых для решения исходной задачи:и то существенного увеличения затрат на расчет траекторий может не произойти.
Здесь следует отметить также, что возникающее в результате преобразования изменение коэффициента сноса в системе СДУ может существенно повлиять на время достижения диффузионным процессом границы области. Поэтому в процессе минимизации время вычисления целевой функции и её градиента может сильно различаться на разных этапах минимизации, так как параметры преобразования меняются при каждом их вычислении. По этой причине в общем случае трудно оценить выгоду от минимизации параметров по сравнению с простым увеличением объема выборки.
Мы приведем сравнение времени вычисления оценки математического ожидания функционала (4.41) для задачи без преобразования и задачи, которая получается в результате преобразования для набора параметров: «о = 1, ot1 = О, (12 = 0, «з = 0. Этот набор параметров соответствует тождественному преобразованию и = и, но при этом для вычисления оценки требуются затраты, связанные с преобразованием. Все расчеты оценок математического ожидания функционала (4.41) проводились с использованием процессора Intel Core 2 Duo Т7100 1.8ГГц.
Расчеты показали, что среднее время вычисления оценки (при объеме выборки 1000) без использования преобразования получилось равным 209 сек., а с использованием преобразования 355 сек. Для вычисления градиента целевой функции численно решалась система СДУ (4.40) и система СДУ, которая получается из нее в результате дифференцирования по параметрам функции V. Среднее время вычисления градиента составило 853 сек.
Если не использовать преобразование, то для определения оценки с такой же точностью, с какой была получена оценка в результате минимизации, надо увеличить объем выборки до 39000 траекторий. На эти вычисления потребовалось 8373 сек. При этом была получена более точная оценка только в одной точке х 1 . Для вычисления оценки с минимизацией параметров потребовалось 6388 сек. Как показали расчеты, определение параметров в результате минимизации относительно этой точки позволило определить более точные оценки в других точках области, значительно удаленных от нее, причем эти оценки были полу 134 чены при объеме выборки 1000 траекторий.
В данном разделе изложен вариант метода минимизации дисперсии функционала диффузионного процесса для случая отражающей границы области. Предполагается , что задана выпуклая область G С R с границей dG, на которой задано условие отражения в направлении внутренней нормали.
Также в данном разделе получена формула для предельного (при убывании шага в методе Эйлера) значения дисперсии требуемой оценки, аналогичная формуле для дисперсии, данной в разделе 4.1.1. Для уменьшения дисперсии вводится параметризованное преобразование параболической краевой задачи с последующей минимизацией выборочной дисперсии оценки. В результате минимизации получается набор параметров, которому соответствует функционал такой, что при вычислении его математического ожидания получается оценка с гораздо меньшей дисперсией.
Пусть х Є G:t Є (0,Т). Предположим, что задан диффузионный процесс X, который описывается системой СДУ S S S Xs = х + / a(v,Xv)dv + / a(v,Xv)dWv + / n(Xv)dKv, (4.47) t t t (4.48) s Ks= f XdG(Xv)dKv, (4.49) t 135 где п - единичный вектор внутренней нормали; Кш - локальное время процесса X на dG. Пусть с, ту, /, 7 : QT - R некоторые достаточно гладкие функции. Определим функции У, Zm следующим образом S S У, = ехр( f c(v,Xv)dv+ frj(v,Xv)dKvy (4.50) t t s s Zs = f f(v,Xv)Yvdv+ [ y(v,Xv)YvdKv. (4.51) t t Мы рассматриваем оценку математического ожидания следующего функционала 4)(XT)YT + ZT (4.52) Известно, что при (t, х) Є QT математическое ожидание этого функционала u(t,x) = Et (XT)YT + ZT ) (4.53) совпадает со значением решения в точке (t, х) следующей краевой задачи для параболического уравнения с заданным условием при t = Т Lu + f(t,x) = 0, (t,x)eQT, (4.54) и(Т,х) = (р(х), х Є G, (4.55) ди — + r](t}x)u + y(t}x) = 0, (t,x)eST, (4.56) \J I v d d где Lu = f + 12 J2 hj(t, x)JL- + J2 ai(t, ж)# + c(t, x)u\ матрица В с элемен 1,3=1 г=1 тами bij удовлетворяет условию В = а о . Заметим, что при введении замены переменной времени t = T—t задача (4.54)—(4.56) становится начально-краевой. Обозначим i) = LP(XT)YT + ZT- Также, как для случая поглощающей границы, в данном разделе получена формула для дисперсии #, из которой следует, что значение дисперсии # зависит от ( ) и У в точках (v,Xv): v Є [0,Т]. Уменьшение значений этих случайных величин в среднем на траекториях Хш приводит к уменьшению дисперсии і). Для уменьшения дисперсии оценки математического ожидания (4.53)предлагается использовать преобразование краевой задачи (4.54)—(4.56) в результате которого первоначальное решение и и новое решение и связаны сотношением u = uV , (4.57) где V(x) 0. Затем рассматривается краевая задача относительной н решается соответствующая ей система СДУ, с которой связывается случайная величина $ = ф(Хт)Ут + %т. Понятно, что при таком преобразовании краевой задачи происходит изменение оценки функционала и ее дисперсии. Поэтому желательно, чтобы в результате такой замены выполнялось неравенство \R$\ \R&\. После того, как получена оценки для й, оценка для исходной задачи определяется по формуле (4.57).
Функция V выбирается, зависящей от некоторого набора параметров, поэтому решение преобразованной краевой задачи тоже зависит от этих параметров. Поскольку в формуле для дисперсии оценки присутствует градиент решения краевой задачи, то и дисперсия оценки решения преобразованной краевой задачи также зависит от этих параметров. Затем с использованием стохастического квазиградиентного алгоритма с переменной метрикой [61] определяется набор параметров, которому соответствует минимальное значение квадрата коэффициента вариации #. Градиент целевой функции вычисляется на основе решения системы СДУ (4.47), (4.49) и системы СДУ, которая получается из нее в результате дифференцирования по параметрам. Для численного решения СДУ используется метод Эйлера.
Вариант метода для условия отражения на границе
Был проведен расчет с использованием СДУ значений температуры в точке, находящейся в центре шестиугольника и на одинаковом расстоянии от верхней и нижней границ панели, на интервале времени от 0 до 250с. При этом панель рассматривалась как неоднородный материал. Для сравнения был проведен расчет значений температуры в той же самой точке, но при этом пластина считалась однородной с усредненным значением температуропроводности, полученным на основе инженерных расчетов. Полученные результаты этих расчетов приведены в виде графиков на рис. 3.
Отметим, что в программу расчета была заложена декартова система координат с началом в центре границы нижней стороны панели, ось Z направлена перпендикулярно панели, а оси X и Y - параллельно панели. Границы шестиугольных сот образованы сетью параллельных оси X прямых и прямых под углом 7г/3 к оси X на расстоянии друг от друга, которое определяется размером сот. В ходе расчетов благодаря циклической структуре расположения прямых, образующих соты, на каждом шаге определяется расстояние до ближайшей прямой. Таким образом определяется ориентация траектории относительно расположения сот.
В процессе вычислений использовалось функция сглаживания7d( -2 ), где 7d - нормирующий коэффициент. Величина шага в методе Эйлера бралась равной 10-6 с таким расчетом, чтобы вероятность попасть за один шаг из одной соты в другую была мала; величина р была задана равной 10-5. Моделировалось 4000 траекторий случайного процесса. При этом объеме выборки ширина 95% доверительного интервала не превышала 1% от значений температуры во всех точках интервала времени, для которого проводился расчет.
Из графика на рис. 3 видно хорошее совпадение значений температур, полученных в результате вычислений с использованием СДУ с результатами расчета соответствующих значений как для однородного материала.
В работе [92] для такого же образца сотовой панели что и в предыдущем примере приведены данные расчета температуры с граничными условиями третьего рода. Для расчета были использованы входные данные, полученные на основе измерений физических данных в течение первых 150с полета самолета в холодном климате. На границе было задано условие hj 3 = 0 , (5.46) xdG где out(t)i если х находится на верхней границе ain(t), если х находится на внутренней границе —otout(t)uout(t), если x находится на верхней границе, ain(t)uin(t), если х находится на внутренней границе. с out, ain - соответственно внешний и внутренний коэффициенты теплообмена, которые рассчитываются на основе экспериментальных данных по методикам, описанным в работах [14], [32]; uout: щп - функции температуры окружающей среды соответственно с наружной и внутренней сторон обшивки.
Таким образом, расчет температуры в заданной точке сотовой панели осуществлялся путем оценки решения в этой точке краевой задачи (5.1), (5.3), (5.46) с использованием численного решения системы СДУ с нулевым вектором сноса Хг = х+ I (7{т){Т -v,Xv)dWv+ f nB(T-v,Xv)d\kv\, (5.47) где пв - единичный вектор внутренней конормали (пв = ш т); 1 1 = I х(Хг dG)d\kr\ - неотрицательный случайный процесс, который возрастает только тогда, когда процесс X. находится на границе. Оценка температуры внутри пластины в точке х в момент времени s определялась по формуле (р{Хт) exp ( rj(T-v, Xv)dv) + / 7(T - v, Xv)dv . (5.48 Для расчёта теплового состояния панели была разработана параллельная программа на языке Fortran 90. Распараллеливание в программе осуществляется по схеме ведущий-ведомые (Master-Slave). В этой схеме одно вычислительное ядро считается главным, и оно распределяет весь объём работы по моделированию случайных траекторий по всем ядрам, участвующим в работе. По окончании моделирования всех траекторий все ядра передают ведущему ядру полученные результаты расчетов для вычисления математического ожидания функционала, дающего оценку температуры. При написании параллельной программы использовалось программное обеспечение Intel MPI, Version 4.1. Моделирование траекторий случайного процесса осуществлялось с использованием параллельного датчика гауссовских случайных величин из библиотеки Intel MKL [15]. Вычисления проводились в Сибирском Суперкомпьютерном центре на гибридном кластере HKC-30T+GPU с использованием 36 4-х ядерных процессоров Е5540 2,53 GHz. Расчеты температуры проводились вблизи нижнего края панели в точке (жо, yoi ZQ) С координатами XQ = уо = 0, ZQ = 0,0001м. Величина шага в методе Эйлера бралась равной 5 10 7 и объемом выборки 4000 случайных траекторий.
При этом величина доверительного интервала температуры сотовой конструкции не превышала 1, 2 в нестационарных и 0,7 К в стационарных условиях при пературы сотовой конструкции с помощью предложенного в данной работе метода. Затем аналогичные расчёты были выполнены для однородной (гомогенной) структуры с экспериментально полученными значениями эффективной теплопроводности и температуропроводности. Полученные результаты рис. 4 показывают, что значения температуры для сотовой конструкции и однородной структуры находятся в хорошем согласии, когда процесс передачи тепла близок к стационарному.