Адаптивно-статистические методы в некоторых задачах вычислительной механики

Отметим также, что функционал качества работы системы определяется не только управляющим воздействием р (х), но и интегрируемой функцией /(х), вносящей фактор неопределенности в процесс управления. Поэтому решение оптимизационной задачи по выбору рорг (х) не приводит к физически реализуемому управлению, т.к. ?>ор*(х) = рор1,{-х., Тем не менее, получение явного выражения для /} может быть использовано для организации процесса адаптивной подстройки вычислений.

При условиях (0.5) оптимальная в смысле минимума функционала (0.6) плотность р*(х) удовлетворяет следующему нелинейному интегральному уравнению: х) ^ + Л?(х) + Ш у)-Ч)2' (0″ 7) где Л выбирается из условия нормировки плотности р (х), а /л 6 Ят-1-из условий:

I (р (х)р* (х)с1х = 0. о

Полученное уравнение для р*(х) дает возможность построения стратегии управления, необходимой для синтеза физически реализуемых адаптивных управлений.

Физическая нереализуемость плотности р*(х) состоит в том, что она существенным образом зависит от функции /(х) и параметров в*, также определяемых через значения /(х). Стратегия адаптивного управления процессом вычисления будет состоять в том, что на основе предыдущих к серий вычислений интегрируемой функции /(х) в отдельных точках строится ее аппроксимация Д (х) во всей области И и дается оценка 9^ и значениям параметров 9 и 7 в соответствии с алгоритмом (0.3). Плотность генерирования узлов сетки в следующей к + 1-ой серии вычислений выбирается в соответствии с (0.7) по формуле:

К+1(х) =, 1 /<,(х) 1 ^ - (0.8)

1 у/Х + /¿-Мх) + К '

Различие алгоритмов адаптивного управления в изложенной постановке определяется различием способов аппроксимации функции /(х). В диссертации рассматриваются два способа аппроксимации /(х) в кусочно-постоянно виде и на основе регрессионной модели (0.2).

Третья глава посвящена адаптивной оптимизации квадратурных формул для вычисления одномерных интегралов на случайной сетке. Результаты, полученные в этой главе, применяются в главе 4 при решении задачи об адаптивном выборе одномерной сетки в методе конечных элементов.

Пусть требуется вычислить одномерный интеграл: ъ

•7 = / ?{х)йх. (0.9) а

Для приближенной оценки интеграла (0.9) используются случайные квадратурные формулы общего вида: n га

JN=Л (х (к+1) — Х (к))? (0.10)

Аг=0 г=1

Здесь ?(?) — случайные узлы интегрирования, причем, а — х^ < х^) < х{к+1) < х{ы+1) — Ъ, к — 1,., ЛГ — вариационный ряд, составленный из независимых реализаций х, случайной величины распределенной с некоторой плотностью р{х) такой, что: р{х) >0, Iр (х)йх = 1, (0.11) некоторые числа, такие что А,' > 0л = = 1, € 1

1х{к), х{к+1)] •

Рассматриваются следующие частные случаи квадратурной формулы (0.10):

1) Формулы 'левых', 'правых' прямоугольников:

А =? — Х (к)), (0.12) к=0

А =? 1(х (к+1)){х (к+1) — Х (к)) — (0.13) к-о

2) Формула трапеций:

М 1

А =? 2(/(х{к)) + 1(х[к+1)Жх (к+1) — х (к)) — (0.14) к=0

3) Формула центральных прямоугольников: n ,

Е/

Х (к) + х (к+1) 2 (^(А+1) — х (к))

0.15)

Аг=0

В литературе квадратурные формулы типа (0.10) рассматриваются для случаев детерминированных узлов. При равномерной расстановке узлов известны соответствующие оценки погрешности вычисления интегралов по этим формулам. Например, у Никольского С. М. приводится следующая оценка: если формула (0.10) точна для многочленов степени г — 1 и /(х) 6 ?(г)(М, а, Ь), т. е. функция /(ж) имеет ограниченную константой М на [а, 6] кусочно-непрерывную производную порядка г, то — 3 <

Обычно наилучшую квадратурную формулу выбирают на определенном классе функций Р, минимизируя остаточный член квадратурной формулы #(/). Узлы квадратурной формулы выбираются из условия минимума величины Бир|Я (/)|. Так как сетка должна быть наилучшей для всех / € Р, то чаще всего оказывается, что узлы х^ -равномерны.

Эта оценка является сильно завышенной на достаточно широком подклассе Р' класса Р (Р1 < Р). Поэтому предлагается вычислять интеграл от функции /, минимизируя оценку погрешности этого вычисления, учитывающую специфику /(х), т. е. адаптивно по отношению к /(#). Для этого предлагается построить оценки остаточного члена Ы^), которые учитывали бы специфику интегрируемой функции / и удовлетворяли неравенству:

Минимизировать оценку погрешности вычисления интеграла Гдг (/) в детерминированном виде практически сложно. Однако, если узлы случайны, то можно получить функционал — М{Т^}, зависящий от функции плотности распределения случайных узлов р (х). Решив вариационную задачу минимизации этого функционала, получаем оптимальную плотность распределения узлов сетки р*(х, /), которая не может быть непосредственно реализована из-за зависимости от /(х).

Ь-а)г+1сгМ

Ыг где сг > 0. вР J-JN < Тдг (/) < 8ИР |Д ^ е^

Основной идеей работы является синтез квадратурных формул на случайной сетке адаптивно приспособленной к конкретной интегрируемой функции /.

Для функций /(х)? С[а] получены оценки Тдг остаточного члена квадратурных формул (0.12), (0.13) (г = 1) и формул (0.15), (0.14) (г = 2). Оценки Тм представляют собой суммы следующего вида: n тм = Е (0.16) 0 где узлы = 0,.,./У" + 1 такие, что Ж (0) — — Ь и х^ < £(г-+1) и функции #(?(г),?(г+1)) для каждой из квадратурных формул (0.12)-(0.14) определяются соответственно выражениями:

1)л (®-(,-),^+1))= / (*(.-+!)-0(0.17)

Щх{{), х{{+1))= I (0.18) *(0

2)^(г-), х (г-+1)) = - / (0.19)

3)Я (*<0,*(Ш)) = /+1)|/" (01 + / |Г (01 ^^, 20) где г, — = «<'>+'<'+ч.

Вывод выражения для М{Тдг} представляет собой самостоятельную задачу, поскольку узлы Хф, входящие в (0.10) не являются стохастически независимыми величинами и вывод выражения для М{Т^} требует привлечения методов теории порядковых статистик.

Сумма Га, — является частным случаем суммы 5дг общего вида: n вм = (0.21)

1=0 где узлы х^, г — 0,.,^ + 1 такие, что = = Ь и аг^ < Ж (г-+1),

Я (ж (г), ж^г+1))-некоторая функция.

Имеет место следующая теорема! :

Теорема1 Пусть H (x^, x^i+iy)-ограниченная на [а, б]2 функция, тогда математическое ожидание Sn определяется выражением: ь ъ

M{SN} = Nj Н (а, х) [1 — P (x)]N~1 p (x)dx + N j H{x, b) P (x)N~lp (x)dx + а, а b b

N (N -1) ff H (x, y) [P (x) + 1 — P (y)]N-2p (y)p (x)dydx1 (0.22) a x x f, а где P (x) = fp (?)d?- функция распределения случайной величины х, М{.} -символ математического ожидания.

Пусть функция Н (х, у) такая, что

1)Н (х, у) — кусочно-дифференцируемая функция своих аргументов,

2)Н (х, х) =0. (0.23)

Тогда для M{Sn} справедливо также представление: ь ь

M{SN} = SI f Щх (х, у)[Р (х) + 1 — P (y)f dydx, (0.24) а x b где S = jH’y (x, y) x—y dy

С помощью теоремы 1 получено аналитическое представление функционала AIn = M{TN}.

Следствие.Математическое ожидание TV определяется выражением: ь ь

AIn = - f f [P (x) + 1 ~ P (y)fdydx. (0.25) а x

Ниже приводятся функционалы Д/У для каждой из квадратурных формул (0.12)-(0.14).

1) Для формулы «левых», «правых» прямоугольников: ь ъ д4? = / I/'Wl / [Р (*) + 1 — P (y)f dydx, (0.26) а х b b

AIP =II If'(y) [P (x) + 1 — P (y)]N dxdy. (0.27) a x

2) Для формулы трапеций: b Ь i 2

А/^ = j |/" (01 // + 1 «dxdydt. (0.28) i а

3) Для формулы центральных прямоугольников:

Ъ Ъ ?

Р (аг) + 1 — dxdyd?. (0.29)

В работе исследовалась проблема оптимизации полученных функционалов по функции плотности р (х) или по функции распределения Р (х). Решить эту задачу аналитически в общем случае не удается, ввиду сложности минимизируемых функционалов и сложных ограничений, накладываемых на функцию Р (х). Однако можно получить выражение для плотности, минимизирующей функционалы (0.26)-(0.29) при числе случайных узлов N" 1.

Справедлива следующая лемма:

Лемма1 Пусть выполнены следующие условия:

1) (г — 1) — максимальная степень многочленов, для которых точна случайная квадратурная формула (0.10);

2) f (x) е W$b] для формул (0.12)-(0Л5) и f (x) G wffi для формулы

0.Ц);

3)плотностъ р (х) > 0- за исключением, быть может, тех точек, в которых f^(x) = 0 и р (х) — кусочно-дифференцируемая функция на [а, Ъ].

Тогда сг brf{r)(x) I (1 1 где сг-некоторый числовой коэффициент, не зависящий от N и f.

Оптимизация функционала (0.30) при числе случайных узлов N «1 по функции плотности р (х) дает следующее выражение для асимптотической оптимальной плотности р (х): р (х) = ср^ |/М (*)|, (0.31) где ср = ь 1 -нормировочный коэффициент. а

Непосредственно воспользоваться плотностью р{х) для расстановки узлов нельзя, поскольку неизвестным является значение нормировочного множителя ср, вычисление которого представляет собой по сложности такую же задачу, что и вычисление самого интеграла 3. Однако, зная аналитическое выражение для оптимальной плотности, можно построить адаптивную процедуру оптимизации случайной сетки.

С другой стороны, функционал (0.25) может быть оптимизирован численно на классе кусочно-постоянных плотностей. При этом не требуется точного учета вклада функции Щ (х, у) в функционал, поскольку достаточно знать только тенденции концентрации случайных узлов.

Для формул трапеций и прямоугольников определена также асимптотическая оптимальная плотность из класса кусочно-постоянных функций.

Доказана теорема о скорости сходимости случайных квадратурных формул (0.10).

Теорема2. Пусть случайная квадратурная формула (0.10) точна для многочленов степени г — 1 и функция /(х)? И^Г)(МГ, а, 6). Пусть плотность р (х) такая, что р{х) > тр > 0, тогда справедлива следующая оценка:

М{АЗ} <-^-, с < оо. (0.32) т/ П (^ + г) г=1

Для формул трапеций и прямоугольников показано, что при числе случайных узлов N >> 1 узлы х*, доставляющие минимум детерминированным оценкам XV и оптимальная асимптотическая плотность р (х) связаны соотношением: х^М{х[{]/р{х)=р (х)}. (0.33)

В третьей главе также приводятся результаты численных экспериментов, в ходе которых методами (0.12)-(0.14) вычислен ряд тестовых интегралов.

Результаты вычисления показали эффективность использования оптимальной плотности (0.31) по сравнению с применением равномерной плотности, а также по сравнению с детерминированной равномерной расстановкой узлов интегрирования. Посчитаны величины характеризующие относительное отклонение усредненной погрешности вычисления интеграла по случайной квадратурной формуле с оптимальным распределением узлов по сравнению с погрешностью вычисления интеграла по детерминированной квадратурной формуле с равноотстоящими узлами. В частности, при вычислении интеграла 3 — / еГх][Х^х методом трапеций с адаптацией плотноо сти величина ЪЗ^ъ для X = 7 составила приблизительно 23%, для, А = 10 -приблизительно 66%. В случае, когда в качестве узлов брались величины, определенные формулой (0.33) ?/50 ~ 82% и 92% соответственно.

Четвертая глава посвящена адаптивной оптимизации случайной одномерной сетки метода конечных элементов при решении задач теории упругости, термоупругости.

При решении практически важных инженерных задач методом конечных элементов одной из важных проблем и, по существу, первым шагом по пути к решению задачи является дискретизация области, поскольку плохое или несовершенное разбиение будет приводить к ошибочным результатам, даже если остальные этапы метода осуществляются с достаточной точностью.

Дискретизация области (тела) включает задание числа, размеров и формы подобластей, которые используются для построения дискретной модели реального тела. При этом возникает двоякая ситуация. С одной стороны, элементы должны быть выбраны достаточно мелкими, а с другой стороны, применение достаточно крупных элементов сокращает вычислительную работу. Нужно иметь некоторые общие соображения об окончательных значениях, с тем, чтобы можно было уменьшить размеры в тех областях, где ожидаемый результат может очень сильно меняться (большие величины градиентов) и увеличить их там, где ожидаемый результат почти постоянен.

В литературе (например, Сегерлинд) указывается, что искусство разбиения области на подобласти зависит от имеющихся инженерных навыков и не имеет теоретического обоснования.

В данной работе выбор сетки определяется следующими соображениями:

Пусть ПЩж)} - потенциальная энергия системы. Тогда если и*(х) -решение, то

П{м*(ж)} < П{и (ж)}, (0.34) где и (х) — любая функция, удовлетворяющая заданным граничным условиям. В МКЭ аппроксимация решения й (х) представляет собой кусочно-полиномиальную функцию, такую что на каждом элементе аппроксимирующий полином имеет следующий вид :

M (ar) = u^(e)W* (0−35)

Здесь т — число узлов элемента, зависит от вида используемого элемента, е — номер элемента, ф^ -вектор-функция, компонентами которой являются координатные функции узлов элемента, um — вектор узловых значений искомой величины, Т обозначает операцию транспонирования.

Узловые значения щ,., Идг, где N — общее число узлов конечно-элементной модели, должны быть теперь «отрегулированы» таким образом, чтобы обеспечивалось «наилучшее» приближение к истинному распределению перемещения. В МКЭ это регулирование осуществляется, применительно к задачам теории упругости, путем минимизации потенциальной энергии по неизвестным компонентам вектора идг после подстановки в ее выражение аппроксимации й (х), т. е. идг = arg min П (и/г, хдг)? (0.36) где хд? = [ж (1), ., Ж (^)]г-вектор упорядоченных узловых значений координат.

Процедура минимизации сводится к решению систем линейных алгебраических уравнений относительно узловых значений щ, В случае если заданы числовые значения координат узлов х^,., #(jv), это позволит сразу получить решение.

В данной работе предлагается координаты узлов не задавать заранее, а рассматривать их реализациями случайной величины ж, распределенной с некоторой плотностью вероятности р (х) такой, что

Ip (x)dx = 1, р (х) >0, х? D (0.37) D и оптимальную плотность распределения находить совместно с узловыми значениями перемещения щ,., и^ из условия:

Рор*(ж), и5г) = arg min М{П (и^, хдг)}, (0.38) p{x), MN где М{-} - символ математического ожидания.

Процедура определения оптимальной плотности может быть осуществлена различными способами.

В ходе одного из них узлы, входящие в выражение (0.36), рассматриваются как параметры. Тогда узловые значения u,., u*N, найденные из условия (0.36), будут зависеть от координат узлов, т. е. vl*n = U*n (Xn)> (0.39)

В этом случае оптимальная плотность popt{x) подлежит определению из условия:

PoPt (x) = argminM{n (uJr (xjv), Xtf)}. (0.40) р{х)

С плотностью popt (x) разыгрываются узлы x*N и в качестве искомых узловых значений перемещения берутся величины и^(х^).

В тех случаях, когда форма конечного элемента учитывается не точно, т. е. аппроксимируется каким-либо образом, условие (0.34), вообще говоря, не выполняется и оптимальная плотность должна быть определена из условия: р^{х) = aigmmM{U (u*N{xN), xN)} (0.41) р (х) где П* = П{и*}.

В диссертации рассматриваемый метод иллюстрируется решением задачи о нахождении поля перемещений в стержне переменной площади поперечного сечения F (x), подверженного осевой нагрузке Q.

Задача решалась симплекс-методом. Рассматривались различные конечно элементные модели балки в зависимости от способа аппроксимации Fi (x) площади поперечного сечения стержня на элементе [?(г),?(г-+1)]- Так, для моделей М/, Mt соответстветственно Fi (x) = F (x^) и F{(x) = ,

Рассмотрен также случай точного учета формы стержня (модель М0).

Для каждой из вышеуказанных моделей доказано, что M{U (VLn (xn), Xn)} = П* + AOfe-F} • Оптимальная плотность определялась из условия:

Popt (x) = arg min |ДП{р, F}.

Определена асимптотическая оптимальная плотность р (х).

Численный эксперимент проведен для стержня длины L = 1 с площадью поперечного сечения F (x) = -F (0)(1 — ?x), где параметр? = 1 — yj^j «

Посчитана величина средней погрешности |Дг/(ж)| = — й (х) оценки перемещения и (х) при равномерном и оптимальном распределении случайных узлов. Оказалось, что применение оптимальной сетки в модели Mi приводит к уменьшению погрешности оценки перемещения и (х) в области х G [0.9,1], т. е. к более точному моделированию перемещения на конце стержня. Этот результат очевиден, поскольку в данном примере оценка для потенциальной энергии и перемещения на конце стержня совпадают с точностью до мультипликативного множителя и минимизация потенциальной энергии автоматически ведет к оптимизации погрешности оценки перемещения на конце стержня.

При точном учете формы стержня (модель М0) величина погрешности перемещения |Дм (я)|, построенного на оптимальной сетке, приблизительно одинакова во всей области стержня, в то время как при использовании равномерной сетки величина погрешности оценки перемещения нарастает к концу стержня и в несколько раз превосходит погрешность оптимального решения.

Получены значения Л^ (индекс указывает на модель стержня) числа детерминированных равноотстоящих узлов, необходимых для достижения той же точности по значению потенциальной энергии, что и при N оптимальных узлах. Так, для достижения той же точности по значению оценки потенциальной энергии для модели М/ требуется число равноотстоящих узлов IVI «1.5N при (3 — 0.9 и iг^ «4.5ЛГ при ?3 = 0.99, для модели М1 -ЛГ* «1Ш при [3 = 0.9 и Щ «7ДГ при ?3 = 0.99.

Проведен сравнительный анализ стохастических и детерминированных решений, построенных по детерминированной равномерной сетке. Выяснилось, что моделирование с использованием равномерной плотности распределения случайных узлов в среднем приводит к более высокому значению погрешности оценки энергии по сравнению с погрешностью, получающейся при детерминированной расстановке узлов равномерным образом. То же самое можно сказать и о моделировании с использованием оптимальной плотности. При моделировании с оптимальной плотностью наблюдается эффект уменьшения величины зтой погрешности по сравнению с детерминированной равномерной расстановкой узлов при значениях параметра ?3 > 0.95.

Предложены различные алгоритмы усреднения стохастического решения й (х), позволяющие уменьшить норму величины 'Аи (х).

IS

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

1. Разработан адаптивно-регрессионный подход к вычислению многомерных интегралов методам Монте-Карло. Подынтегральная функция аппроксимируется функцией регрессии, параметры которой определяются методом наименьших квадратов (МНК).

2. Получено аналитическое выражение функционала дисперсии оценки интеграла регрессионным методом. Показано, что величина дисперсии может быть значительно уменьшена по сравнению с традиционным методом Монте-Карло за счет подходящего выбора базисных функций.

3. Получено с точностью до параметров аналитическое выражение для плотности минимизирующей функционал дисперсии.

4. Предложены 2 алгоритма адаптации плотности р{х, /}.

5. Разработаны алгоритмы адаптивной оптимизации вычисления одномерных интегралов по квадратурным формулам со случайными узлами.

6. Получены функционалы оценок погрешностей вычисления интегралов. Эти функционалы зависят от плотности распределения случайных узлов и от характеристик интегрируемой функции.

7. Получено аналитическое выражение для асимптотической оптимальной плотности, минимизирующей функционал погрешности. Эта плотность зависит от характеристик интегрируемой функции. Построена ада

О и < птивная процедура оптимизации погрешностей вычислении интегралов, эф-фективнотсь использования которой подтверждена численными экспериментами.

8. Предложен подход к рациональному выбору случайной сетки метода конечных элементов при рассчете одномерных стержневых систем. Этот подход заключается в том, что плотность распределения случайных узлов определяется из условия минимума математического ожидания потенциальной энергии, выраженной через аппроксимации искомого решения.

9. Получены основные соотношения метода для стержня переменной площади поперечного сечения, подверженного осевой нагрузке. Численный эксперимент проведен на примере стержня с линейно меняющейся площадью поперечного сечения. Результаты численного моделирования показали достаточно высокую эффективность использования адаптивно-оптимальной случайной сетки МКЭ.

Публикации по теме диссертации

1. Кульчицкий О. Ю., Иванов В. М., Бутенина Д. В. Вычисление интегралов методами статистического моделирования с адаптацией. Труды СПбГТУ, сер. Механика и процессы управления, N458, Изд. СПбГТУ, С.-Петербург, 1995, с.151−161.

2.Butenina D.V., Ivanov V.M., Kul’chitsky O.Yu. The adaptive optimization of nodes distribution in multivariate quadrature formulas. International Conferena «Optimization of Finite Element Approximations», St. Petersburg, Russia, 1995, p.46−47.

3. Бутенина Д. В. «Оптимизация сетки метода конечных элементов в задаче о растяжении балки переменного сечения.» 5-я Международная конференция женщин-математиков «Математика. Экономика.» Тезисы докладов. 26 июня — 1 июня 1997. Ростов на Дону, с. 7.

4. Бутенина Д. В., Кульчицкий О. Ю., Иванов В. М. «Оптимизация случайной сетки метода конечных элементов в задаче о растяжении стержня переменного сечения.» Международный семинар. Нелинейное моделирование и управление. Тезисы докладов. Самара, 24−27 июня 1997 г., с. 158.

5. Бутенина Д. В., Кульчицкий О. Ю. «Адаптивно-статистический подход к решению задач вычислительной механики.» Труды IV Международной конференции женщин-математиков «Математика. Моделирование. Экология» Волгоград, 27−31 мая 1996 г., т. 4, вып. 1, Изд.Изв.Вузов, радиофизика, Нижний Новгород, 1997, с. 65−72.

6. Иванов В. М., Кульчицкий О. Ю., Бутенина Д. В. Адаптивное управление вычислительными процессами. Труды СПбГТУ, сер. Механика и процессы управления, N467, Изд. СПбГТУ, С.-Петербург, 1997, с.54−60.

Заключение