Диссертация (Бессеточные методы Монте-Карло решения краевых задач для уравнений в частных производных), страница 9

< tm = t,и Yt0 , Yt1 , . . . , Ytm — эйлерова аппроксимация диффузионного процесса, py (Y ) —значение плотности конечномерного распределения этого процесса на построенной реализации Y , наконец, ζ — несмещенная оценка плотности конечномерного распределения диффузионного процесса, тогда выражение ζg(Ytm )/py (Y ),очевидно, будет несмещенной оценкой для функционала Eg(Xt ).

В работе построены оценки ζ, имеющие конечные моменты при достаточно малом шаге вметоде Эйлера. Результаты работы [76] также представлены в [12]. В работах[77,79,80] рассмотрены алгоритмы метода Монте-Карло для вычисления другихфункционалов на траекториях диффузионных процессов.Наконец, для случая дифференцируемых коэффициентов уравнения,можно получить, используя формулы Грина, интегральное уравнение дляu(x, t) и решить его методом Монте-Карло. Такие алгоритмы рассмотрены вработах [43] и [44] для уравнений, главной частью которых является операторЛапласа.

На уравнение с переменной матрицей старших коэффициентов онинепосредственно не переносятся.Мы будем строить несмещенные оценки решения задачи Коши и функционалов от него на траекториях случайных блужданий. В основе построенийлежит формула (2.2.2), в которой каждый интеграл оценивается по одному случайному узлу. Фундаментальное решение является функционалом от решенияинтегрального уравнения Вольтерра (2.1.7), поэтому оценка для него находитсяпо схеме Неймана-Улама [11].Для уравнений с гладкими коэффициентами само решение задачи Кошиудовлетворяет интегральному уравнению Вольтерра со слабо полярным ядром,к которому также применима схема Неймана-Улама.622.2.1Несмещенные оценки решения задачи КошиИз представления (2.2.2) для решения задачи Коши и формулы (2.1.6)для фундаментального решения следует, что решение задачи Коши распадаетсяв сумму четырех потенциалов:u(x, t) = u1 (x, t) + u2 (x, t) + u3 (x, t) + u4 (x, t) =Zt=ZZ0 (x − y, y, t, τ )f (y, τ )dy +dτ0RnZdτ0Z0 (x − y, y, t, 0)ϕ(y)dy +RnZt+ZZZ0 (x − z, z, t, λ)Q(z, y, λ, τ )dz  f (y, τ )dy +dλRnZtτZ+RnZtdλRnZ0Z0 (x − z, z, t, λ)Q(z, y, λ, 0)dz  ϕ(y)dy.

(2.2.6)RnВ силу неравенства (2.1.5) для несмещенного оценивания u1 (x, t) иu2 (x, t) можно использовать плотность Z1 . Нормально распределенный случайный вектор Y , имеющий плотность Z1 (x − y, t − τ ) можно моделировать поформуле (2.2.3), которая теперь примет видpY (x, t, τ ) = x + 4µ(t − τ )γω.(2.2.7)ПустьY = Y (x, t, tθ),а функцияZ0,Z1тогда, аналогично (2.2.4), несмещенными оценками u1 (x, t) и u2 (x, t) будутW =ξ1 (x, t) = t · W (x − Y, Y, t, tθ) · f (Y, tθ)(2.2.8)ξ2 (x, t) = W (x − Y, Y, t, 0) · ϕ(Y ),(2.2.9)и63соответственно.

Случайные элементы γ, θ, ω в этих формулах должны бытьнезависимы в совокупности. Если функции f и ϕ ограничены, то дисперсиипостроенных оценок конечны, поскольку функция W также ограничена в силунеравенства (2.1.5).Пусть функции f и ϕ — ограничены. Тогда в выражениях для u3 и u4 вформуле (2.2.6) можно менять порядок интегрирования, поскольку повторныеинтегралы сходятся абсолютно, в силу неравенств (2.1.5) и (2.1.13).Рассмотрим функцииZλZv3 (z, λ) =dτQ(z, y, λ, τ )f (y, τ )dy,(2.2.10)Rn0Zv4 (z, λ) =Q(z, y, λ, 0)ϕ(y)dy.(2.2.11)RnТогдаZtZdλu3 (x, t) =0RnZtZu4 (x, t) =dλ0Z0 (x − z, z, t, λ)v3 (z, λ)dz,(2.2.12)Z0 (x − z, z, t, λ)v4 (z, λ)dz.(2.2.13)RnВ силу неравенства (2.1.13), функция v3 (z, λ) — ограничена, аαv4 (z, λ) 6 const ·λ 2 −1 .Умножая уравнение (2.1.7) на f (y, τ ) и интегрируя, получим для v3 интегральное уравнениеZtv3 (x, t) +Zdλ0K(x, z, t, λ)v3 (z, λ)dz +RnZt+Zdτ0K(x, y, t, τ )f (y, τ )dy = 0.

(2.2.14)Rn64Аналогично, подставляя в (2.1.7) τ = 0, умножая его на ϕ(y) и интегрируя поy, получаем для v4 уравнениеZtv4 (x, t) +Zdλ0K(x, z, t, λ)v4 (z, λ)dz +RnZ+K(x, y, t, 0)ϕ(y)dy = 0. (2.2.15)RnИз неравенства (2.1.11) следует, что к уравнениям (2.2.14) и (2.2.15) применимасхема Неймана-Улама [11]. Для ее реализации достаточно выбрать плотностьвероятности перехода обрывающейся цепи Маркова, согласованную с ядромуравнения K(x, z, t, λ).

В качестве такой плотности можно взять, например, α(1 − q)αp (x, t) → (z, λ) =(t − λ) 2 −1 Z2 (x − z, t − λ),α2t 2(2.2.16)где 0 < q < 1 является вероятностью обрыва цепи на текущем шаге, а n2|x − z|2Cexp −C(2.2.17)Z2 (x − z, t − λ) =π(t − λ)t−λпри 0 6 λ < t. При λ > t функция Z2 (x − z, t − λ) = 0. Постоянная C в этихформулах берется из неравенства (2.1.9), которое также влечет согласованностьплотности и ядра интегрального уравнения. Отметим, что вероятность обрывацепи на каждом шаге постоянна, поэтому цепь обрывается с вероятностью 1 исреднее число шагов до обрыва равно q −1 .Для моделирования цепи {(xm , tm )}∞m=1 , стартующей из точки (x0 , t0 ) == (x, t), будем использовать независимые в совокупности изотропные случай∞ных векторы {ωm }∞m=1 и случайные величины {γm }m=1 , имеющие гамма распре-деление с параметром n/2.

Предполагается также, что эти случайные элементынезависимы с выбранными ранее элементами ω, θ, γ.Пусть, кроме того, {βm }∞m=1 последовательность случайных величин,независимых в совокупности с определенными ранее и имеющих на отрезке65[0, 1] бета распределение с параметрами (1, α2 ).

Плотность распределенияp(s) =αα(1 − s) 2 −12(2.2.18)Тогда, при m = 1, 2, . . .tm = tm−1 βm ,pxm = xm−1 + C −1 (tm−1 − tm )γm ωm(2.2.19)Момент обрыва цепи N имеет геометрическое распределение с параметром q,то есть P {N = m} = q(1 − q)m , при m = 0, 1, 2, . . . и независим от траектории.Для построения несмещенных оценок для решения уравненияZtv(x, t) +Zdλ0K(x, z, t, λ)v(z, λ)dz + F (x, t) = 0,(2.2.20)Rnопределим последовательность весовых функций W (0) = 1,W (m) = (−1)m+1 W (m−1)K(xm−1 , xm , tm−1 , tm ),p (xm−1 , tm−1 ) → (xm , tm )(2.2.21)при m = 1, 2, .

. .Стандартными несмещенными оценками ([11]) для функции v(x, t) являются случайные величиныη(x, t) =NXW (m) F (xm , tm )(2.2.22)m=0иζ(x, t) =W (N ) F (xN , tN ).qВ качестве функции F (x, t) в уравнении (2.2.14) беретсяZtF (x, t) =Zdτ0K(x, y, t, τ )f (y, τ )dy,Rn(2.2.23)66поэтому несмещенными оценками для v3 (x, t) будутη3 (x, t) = (1 − q)NXW (m+1) f (xm+1 , tm+1 )m=0и(1 − q)W (N +1) f (xN +1 , tN +1 ).ζ3 (x, t) =qПолагая Y = Y (x, t, tθ), для функции u3 (x, t) получаем две несмещенные оценки:ξ3 (x, t) = t · W (x − Y, Y, t, tθ) · η3 (Y, tθ)(2.2.24)и0ξ3 (x, t) = t · W (x − Y, Y, t, tθ) · ζ3 (Y, tθ)(2.2.25)В уравнении (2.2.15) роль функции F (x, t) выполняет функцияZF1 (x, t) = K(x, y, t, 0)ϕ(y)dy,Rnкоторая неограничена.

Поэтому, аналогичные несмещенные оценки для v4 (x, t)имеют бесконечную дисперсию. Чтобы получить оценки с конечной дисперсией,нужно при оценивании F (x, t) включить особенность ядра K(x, y, t, 0) в плотность . Для этого запишем v4 (x, t) в виде разности v4 (x, t) = v5 (x, t) − F1 (x, t).Тогда для v5 (x, t) справедливо уравнение (2.2.20), в котором функцияZtF (x, t) = F2 (x, t) = −Zdλ0ZK(x, z, t, λ)RnK(z, y, λ, 0)ϕ(y)dydz(2.2.26)Rnуже является ограниченной.

Таким образом,Ztu4 (x, t) = −Zdλ0RnZZ0 (x − z, z, t, λ)K(z, y, λ, 0)ϕ(y)dydz +Rn+ u5 (x, t). (2.2.27)67Для оценки интеграла в (2.2.27) возьмем случайную величину β, имеющую бетараспределение с параметрами ( α2 , 1). ПоложимY = Y (x, t, tβ),y0 = Y,y1 = y0 +pC −1 (tβ)γ1 ω1 ,α K(y0 , y1 , tβ, 0)2W1 = −t β 1− 2,αZ2 (y1 − y0 , tβ)(2.2.28)ξ4 (x, t) = W (x − Y, Y, t, tβ) · W1 · ϕ(y1 ).Случайная величина ξ4 (x, t) — несмещенная оценка интеграла.Функцию u5 (x, t) оцениваем на траекториях цепи Маркова.

Для этого требуется получить несмещенную оценку F2 (xm , tm ). Интегралы (2.2.26) и(2.2.27) имеют одинаковую структуру и одинаково оцениваются. Теперь случайная величина β имеет бета распределение с параметрами ( α2 , α2 ). Построениеоценки завершают формулы:pC −1 (tm − tm β)γm+1 ωm+1 ,py1 = y0 + C −1 tm βγm+2 ωm+2 ,α ααα,β 1− 2 (1 − β)1− 2 ×Wm = −tm B2 2K(xm , y0 , tm , tm β) K(y0 , y1 , tm β, 0)×,Z2 (xm − y0 , tm − tm β) Z2 (y0 − y1 , tm β)y0 = x m +(2.2.29)φm = Wm ϕ(y1 ).Здесь φm — несмещенная оценка F2 (xm , tm ), а α α Z1 αα,= s 2 −1 (1 − s) 2 −1 dsB2 20— бета-функция. По аналогии с оценками для u3 (x, t), получаем несмещенныеоценки для u5 (x, t):ξ5 (x, t) = t · W (x − Y, Y, t, tθ) ·NXW (m) φm ,(2.2.30)m=0W (N ) φNξ5 (x, t) = t · W (x − Y, Y, t, tθ) ·,q0(2.2.31)68которые построены на траекториях цепи Маркова, стартующей из точки (Y, tθ).Почти очевидным следствием включения всех особенностей в плотностивероятностей перехода марковской цепи является конечность дисперсии построенных оценок.

Действительно, дисперсии стандартных оценок в схеме НейманаУлама заведомо конечны ([11]), если ряд Неймана для ядра K 2 /p сходится втом функциональном пространстве, в котором решается интегральное уравне(T )ние (2.2.20). В нашем случае это пространство L∞ (Dn+1 ). Поскольку π n2|K|2c6= c2 ,pα(1 − q) Cто для оператора Вольтерра с ядром K 2 /p выполнено неравенствоn+2−αK 2 (x, y, t, τ )|x − y|2− 6 cc2 (t − τ ) 2 exp −C,t−τp (x, t) → (y, τ )(2.2.32)которое влечет сходимость ряда Неймана.

Отметим, что дисперсии оценок равномерно ограничены при ограниченных f и ϕ.Таким образом, доказанаТеорема 2.2.1 Пусть коэффициенты параболического оператора (2.1.1) приα(T )надлежат классу H α, 2 (Dn+1 ),α < 1, матрица коэффициентов при старшихпроизводных — симметричная, а ее собственные числа лежат в фиксированα(T )ном отрезке [ν, µ] и ν > 0. Пусть функция f (x, t) из класса H α, 2 (Dn+1 ) инепрерывная функция ϕ(x) ограничены. Тогда случайные величины ξ = ξ1 ++ ξ2 + ξ3 + ξ4 + ξ5 и ξ 0 = ξ1 + ξ2 + ξ30 + ξ4 + ξ50 , определяемые формулами (2.2.8),(2.2.9), (2.2.24), (2.2.25), (2.2.28–2.2.31), являются несмещенными оценкамирешения задачи Коши (2.2.1).

Дисперсии оценок равномерно ограничены по переменным x, t.2.2.2Определение постоянных c и CВычислив производные в представлении (2.1.8) ядра K, получим69K(x, y, t, τ ) = −+1[n − T r(A(x, t)A−1 (y, τ ))]Z0 +2(t − τ )1(x − y)0 A−1 (y, τ )[A(y, τ ) − A(x, t)]A−1 (y, τ )(x − y)Z0 −24(t − τ )1−a0 (x, t)A−1 (y, τ )(x − y)Z0 + a0 (x, t)Z0 , (2.2.33)2(t − τ )где функция T r(·) — вычисляет след матрицы, a(x, t) — столбец с компонентамиai (x, t), а операция a0 означает транспонирование матрицы.В качестве C можно взять любую постоянную, удовлетворяющую неравенству4µC < 1.eBe — симметричные матрицы с постоянными элементами, удоПусть A,влетворяющими условиямα|ai,j (x, t) − ai,j (y, τ )| 6 eai,j |x − y|α + ebi,j |t − τ | 2 ,а S — ортогональная матрица, приводящая матрицу A−1 (y, τ ) к диагональномувиду. Тогда,n − T r(A(x, t)A−1 (y, τ )) = T r(S 0 (A(y, τ ) − A(x, t))SS 0 A−1 (y, τ )S) ==nXλj Sj0 (A(y, τ )− A(x, t))Sj =n Xnn XXλj sk,j (ak,m (y, τ ) − ak,m (x, t))sm,j ,j=1 k=1 m=1j=1где Sj — j-й столбец матрицы S, а λj — j-е собственное число матрицы A−1 (y, τ ).Отсюда получаем неравенстваnn1 XX|n − T r(A(x, t)A (y, τ ))| 6|ak,m (y, τ ) − ak,m (x, t))| 6νk=1 m=1nnαα1 XXαe26eak,m |x − y| + bk,m |t − τ | = ec1 |x − y|α + ec2 |t − τ | 2 .