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

Величины ϑ, δ имеют плотности распределения α2 s 2 −1 и1√,2 sсоответственно, а θ распределена равномерно.n+αα 2 Γ(2 )2g1τ)αΓ( n2 )αη1 (y, τ ) = −(T −s !n+αy+2 γ(T − τ )ϑΩ, y, τ ×2!s n+α(T − τ )ϑΩ, τ + (T − τ )ϑ) ,×v y+2 γ2r α 1nη2 (y, τ ) = −(T − τ ) 2 g2 y + 2 γ(T − τ )ϑΩ, y, τ + (T − τ )ϑ, τ ×α2r n×v y+2 γ(T − τ )ϑΩ, τ + (T − τ )ϑ)) ,2s !n+2+αα+1Γ( 2 )α 2n+2+α(T − τ )ϑΩ, y, τ ×η3 (y, τ ) = −(T − τ ) 2h1 y + 2 γnαΓ( 2 )2s !n+2+α×v y+2 γ(T − τ )ϑΩ, τ + (T − τ )ϑ) ,2s !α nn+2(T − τ )ϑΩ, y, τ + (T − τ )ϑ, τ ×η4 (y, τ ) = −(T − τ ) 2 × h2 y + 2 γα2s !n+2×v y+2 γ(T − τ )ϑΩ, τ + (T − τ )ϑ ,22Γ( n+12 )η5 (y, τ ) = −(T − τ )×nΓ( 2 )s !n+1× a0 y + 2 γ(T − τ )δΩ, τ + (T − τ )δ A−1 (y, τ )Ω ×2s !n+1×v y+2 γ(T − τ )δΩ, τ + (T − τ )δ ,2r nη6 (y, τ ) = (T − τ )a0 y + 2 γ(T − τ )θΩ, τ + (T − τ )θ ×2r n×v y+2 γ(T − τ )θΩ, τ + (T − τ )θ .212(2.2.53)80Выбирая с вероятностью16одно из слагаемых в (2.2.53) и умножая его на 6,получаем окончательную несмещенную оценку ζ(y, τ ) для V (y, τ ).Оценки ψm для функционалов Ψm (h) из (2.2.45) будем строить на траекториях неоднородной цепи Маркова {(yk , τk )}∞k=0 .

Начальное состояние цепи(T )выбирается в Dn+1 с некоторой плотностью π(y, τ ), согласованной с функциейf (y, τ ). В оценку ψm при этом записывается дробьf (y0 , τ0 ).π(y0 , τ0 )Далее, последовательно выполняются m переходов, на каждом из которых реализуется оценка ζ(yk , τk ) по независимым в совокупности последовательностям∞∞∞случайных элементов {ϑk }∞k=0 , {δk }k=0 , {θk }k=0 , {Ωk }k=0 . При этом ψm умно-жается на весовой множитель в оценке ζ(yk , τk ), а аргументы функции v определяют следующее состояние цепи. Например, если на шаге k = 1, 2, .

. . , m воценке η(yk−1 , τk−1 ) было выбрано первое слагаемое η1 (yk−1 , τk−1 ), то ψm нужноумножить на2α Γ( n+α2 )− 6(T − τk−1 )×nαΓ( 2 )α2s× g1 yk−1 + 2γk−1!n+α(T − τk−1 )ϑk−1 Ωk−1 , yk−1 , τk−12и перейти в точку (yk , τk ) с координатамиsn+αyk = yk−1 + 2 γk−1(T − τk−1 )ϑk−1 Ωk−1 ,2τk = τk−1 + (T − τk−1 )ϑk−1 .На заключительном этапе, в соответствии с (2.2.45), выбираются время tm , распределенное равномерно на отрезке [τm , T ], и координата xm , распределеннаянормально с плотностью Z0 (x − ym , ym , tm , τm ), а оценка ψm умножается на(T − tm )h(xm , tm ). Применяя формулы (2.2.42–2.2.45), получаем оценки функ-81ционала Φ3 (h):φ3 = −N Xm=11−1−qm−1ψm ,(2.2.54)N −1ψN1.(2.2.55)φ03 = − −1−qqФункционал Φ4 (h) оценивается также по формулам (2.2.54–2.2.55), в которых величины ψm строятся на траекториях марковской цепи, стартующей източки (y0 , 0).

Первоначальное значение оценки ψm равноϕ(y0 ),π(y0 )где y0 распределено в Rn с плотностью π(y), согласованной с функцией ϕ(y).Ограниченность дисперсий построенных оценок вытекает из ограниченности весового множителя в оценке ζ(y, τ ) и доказывается так же, как аналогичное утверждение в “прямой схеме” для первой краевой задачи [52].f 2 (y, τ )— инπ(y, τ )= Rn × (0, T ) . Тогда дисперсии оценок φ3 и φ03 конечны.Теорема 2.2.2 Пусть функция h(x, t) — ограничена, а функция(T )тегрируема в Dn+1Доказательство. Оценим второй момент для φ03 :E(φ03 )2∞X2Eψm=.m−1q(1−q)m=1(2.2.56)Заметим, что2Eψmm−1f 2 (y0 , τ0 ) Y e2=Eζ (yk , τk )(T − tm )2 h2 (xm , tm ),π(y0 , τ0 )k=0e k , τk ) — весовой множитель в оценке ζ(yk , τk ).

Рассмотрим интегральгде ζ(yный оператор, определенный представлением (2.2.51), в котором все интегралы берутся со знаком “+00 , а функции g1 , g2 , h1 , h2 , a0 и выражение a0 (y +p+2 s(λ − τ )Ω, λ)A−1 (y, τ )Ω заменяются на их модули. Пусть K(x, y, t, τ ) — ядро этого оператора, а K m (x, y, t, τ ) — его повторное ядро. Для ядра K(x, y, t, τ )82справедлива оценка (2.1.9) с некоторыми новыми постоянными c1 и C1 , поэтомуряд из повторных ядер∞XB m K m (x, y, t, τ )(2.2.57)m=1e k , τk ) вравномерно сходится для любой постоянной B. Весовой множитель ζ(yоценке ζ(yk , τk ) ограничен в силу ограниченности функций g1 , g2 , h1 , h2 , a0 иpвыражения a0 (y + 2 s(λ − τ )Ω, λ)A−1 (y, τ )Ω. Поэтому, существует постояннаяe k , τk )| 6 (1 − q)B. Тогда,B, такая что |ζ(y2Eψm6ZT0ZTZZTf 2 (y, τ )dτ dydλ dz dt ×π(y, τ )τRnRnλZ× Z0 (x − z, z, t, λ)(1 − q)m B m K m (z, y, λ, τ )h2 (x, t)dx.

(2.2.58)ZRnТеперь, из равенства (2.2.56) следует неравенствоE(φ03 )2 6ZT0ZTZZTf 2 (y, τ )dλ dz dt ×dτ dyπ(y, τ )τRnRnλZ1−q× 2Z0 (x − z, z, t, λ)Q(z, y, λ, τ )T h2 (x, t)dx < ∞, (2.2.59)qZRnгде Q — сумма ряда (2.2.57).Второй момент случайной величины φ3 равен∞X∞∞2XXEψm|ψm |E(φ3 ) 6+2E|ψl |.m−1m−1(1−q)(1−q)m=1m=12(2.2.60)l=m+1Уже доказано, что первое слагаемое в (2.2.60) конечно. Пусть Fm — σ-алгебра,порожденная цепью до момента времени m, а функция h1 (y, τ ) определена равенством:ZTh1 (y, τ ) =ZTdλτZdtλRnZZ0 (x − z, z, t, λ)Q(z, y, λ, τ )|h(x, t)|dx,dzRn83тогдаE∞X!|ψl | | Fml=m+1m−1|f (y0 , τ0 )| Y e|ζ(yk , τk )|h1 (ym , τm ).=π(y0 , τ0 )k=0Теперь, аналогично (2.2.58) и (2.2.59), получаем неравенствоZTZZTZZT∞2f(y,τ)(−1)m ψm X(−1)l ψl | 6 dτ dydλ dz dt ×|Em−1(1−q)π(y,τ)m=1l=m+1τRnRnλZ 0× h1 (z, λ)(1 − q) Z0 (x − z, z, t, λ)Q(z, y, λ, τ )|h(x, t)|dx < ∞.∞XRn2.2.4Задача Коши для уравнений с дифференцируемымикоэффициентамиВ данном параграфе изучаются алгоритмы статистического моделирова-ния решения задачи Коши, полученные путем сведения дифференциальной задачи к интегральному уравнению с помощью формулы Грина.

Для уравнений,у которых коэффициенты при вторых производных постоянные, такие алгоритмы можно найти в монографии [43]. Мы рассматриваем общий случай, когдавсе коэффициенты уравнения являются функциями координат и времени.Если коэффициенты aij (x, t) оператора L имеют ограниченные производные второго порядка, а коэффициенты ai (x, t) — ограниченные производные первого порядка по пространственным переменным, то применив формулу(2.1.35), получим интегральное уравнениеZt Z[v0 (y, τ )f (y, τ ) − u(y, τ )M v0 (y, τ )] dydτ +u(x, t) =0 RnZ+ϕ(y)v0 (y, 0)dy. (2.2.61)Rn84При сделанных предположениях справедливо неравенство (2.1.26). Следовательно, операторZt ZKu(x, t) = −u(y, τ )M v0 (y, τ )dydτ(2.2.62)0 Rnимеет слабую особенность и ряд Неймана для уравнения (2.2.61) сходится, еслифункции f (x, t) и ϕ(x) ограничены.u(x, t) =∞XK i F (x, t),(2.2.63)i=0гдеZt ZF (x, t) = F1 (x, t) + F2 (x, t) =v0 (y, τ )f (y, τ )dydτ +0 RnZ+ϕ(y)v0 (y, 0)dy.

(2.2.64)RnИз формулы (2.1.25) следует, что ядро K1 (x, y, t, τ ) оператора K представимов видеK1 (x, y, t, τ ) = −M v0 (y, τ ) = −M v0 (y, τ ) + M0 v0 (y, τ ) =nnXX∂∂2v0 (y, τ ) −di (y, τ ) v0 (y, τ ) −[−aij (x, t) + aij (y, τ )]=∂yi ∂yj∂yii=1i,j=1[n − T r(A(y, τ )A−1 (x, t))]v0 (y, τ ) +2(t − τ )(y − x)0 A−1 (x, t)[A(y, τ ) − A(x, t)]A−1 (x, t)(y − x)+v0 (y, τ ) +4(t − τ )2d0 (y, τ )A−1 (x, t)(y − x)+v0 (y, τ ) − d0 (y, τ )v0 (y, τ ), (2.2.65)2(t − τ )− d0 (y, τ )v0 (y, τ ) =поэтому для построения несмещенных оценок для функции u(x, t) можно применять как методы §2.2.1, так и методы §2.2.3.85В первом случае, оценки строится на траекториях цепи Маркова с переходной плотностью (2.2.16) и имеют видη(x, t) =NXW (m) F (xm , tm )(2.2.66)m=0иW (N ) F (xN , tN )ζ(x, t) =,q(2.2.67)где N – момент обрыва цепи, имеющий геометрическое распределениеP (N = m) = q(1 − q)m ,m = 0, 1, .

. . .(2.2.68)Весовые коэффициенты вычисляются по формулам W (0) = 1,W (m) = W (m−1)K1 (xm−1 , xm , tm−1 , tm ),p (xm−1 , tm−1 ) → (xm , tm )(2.2.69)при m = 1, 2, . . .. Траектория цепи моделируется по формулам (2.2.19). Постоянная C должна удовлетворять условию 4µC < 1. Окончательные несмещенные оценки для u(x, t) получаются заменой F (x, t) на ее несмещеную оценкуp√tf (x + 2t(1 − θ)Y, tθ) + ϕ(x + 2tY ), где θ распределена равномерно на отрезке [0, 1], а Y — случайный вектор, имеющий нормальное распределение сосредним 0 и матрицей ковариаций A(x, t).Рассмотрим теперь применение методов §2.2.3. Зафиксируем число q(0 < q < 1) и промоделируем геометрическое распределение (2.2.68) с параметром q. Пусть N — полученная реализация, тогдаK N F (x, t)ξ1 (x, t) =q(1 − q)Nиξ2 (x, t) =NXK m F (x, t)m=0(1 − q)mявляются несмещенными оценками для u(x, t).

Для определения несмещеннойоценки для K m F (x, t) последовательно, m раз, выполним процедуру оцени-86вания интеграла (2.2.62), которая аналогична процедуре оценивания интеграла (2.2.46). Формула (2.2.48) теперь примет видZ∞ZtKu(x, t) =dr ×dτ00!2n − T r A(x + rΩ, τ )A (x, t)rexp −rn−1 u(x + rΩ, τ ) +nn24(t − τ )(t − τ )Γ( 2 )(4(t − τ ))ZtZ∞2r2 Ω0 A−1 (x, t)A(x + rΩ, τ )A−1 (x, t)Ω − 1+ dτ drE×n4(t − τ )2 Γ( n2 )(4(t − τ )) 200r2rn−1 u(x + rΩ, τ ) +× exp −4(t − τ )ZtZ∞+ dτ dr ×−1×E00rd (x + rΩ, τ )A (x, t)Ωr2×Eexp −rn−1 u(x + rΩ, τ ) +nn24(t − τ )(t − τ )Γ( 2 )(4(t − τ ))ZtZ∞2rn−1− dλ drEd0 (x + rΩ, τ )u(x + rΩ, τ ) nn ×Γ( 2 )(4(t − τ )) 200r2. (2.2.70)× exp −4(t − τ )0−1Случайный вектор Ω распределен на S10 (x, t) с плотностьюp(x, t, ω) =1p.σn det(A(x, t))|A−1 (x, t)ω|(2.2.71)Аналогами формул (2.2.49) и (2.2.50) будут следующие формулы, в которыхg̃1 , g̃2 , h̃1 , h̃2 являются ограниченными функциями.n − T r A(x + rΩ, τ )A−1 (x, t) = T r [A(x, t) − A(x + rΩ, t)] A−1 (x, t) ++ T r [A(x + rΩ, t) − A(x + rΩ, τ )] A−1 (x, t) =α= g̃1 (x + rΩ, x, t)rα + g̃2 (x + rΩ, x, τ, t)(t − τ ) 2 , (2.2.72)87 0 −1Ω A (x, t)A(x + rΩ, τ )A−1 (x, t)Ω − 1 == Ω0 A−1 (x, t) [A(x + rΩ, t) − A(x, t)] A−1 (x, t)Ω ++ Ω0 A−1 (x, t) [A(x + rΩ, τ ) − A(x + rΩ, t)] A−1 (x, t)Ω =α= h̃1 (x + rΩ, x, t)rα + h̃2 (x + rΩ, x, τ, t)(t − τ ) 2 , (2.2.73)Подставляя эти выражения в (2.2.70) и выполняя замену переменнойs = r2 /4(t − τ ), получаем для Ku(x, t) представление.Ztdτ (t − τ )Ku(x, t) =α2 −1Z∞002α n+α −1exp(−s) ×dss 22Γ( n2 )pp× E g̃1 (x + 2 s(t − τ )Ω, x, t)u(x + 2 s(t − τ )Ω, τ ) +Ztdτ (t − τ )+Z∞α2 −1ds00n12 −1 exp(−s) ×sn2Γ( 2 )pp× E g̃2 (x + 2 s(t − τ )Ω, x, τ, t)u(x + 2 s(t − τ )Ω, τ ) +Ztdτ (t − τ )+α2 −10Z∞02α n+2+α −1ds n s 2exp(−s) ×Γ( 2 )pp× E h̃1 (x + 2 s(t − τ )Ω, x, t)u(x + 2 s(t − τ )Ω, τ ) +Ztdτ (t − τ )+α2 −1Z∞ds001 n+2 −1s 2exp(−s) ×Γ( n2 )pp× E h̃2 (x + 2 s(t − τ )Ω, x, τ, t)u(x + 2 s(t − τ )Ω, τ ) +Zt− 12Z∞dτ (t − τ )+0ds01 n+1 −1s 2exp(−s) ×Γ( n2 )pp0−1× E d (x + 2 s(t − τ )Ω, τ )A (x, t)Ωu(x + 2 s(t − τ )Ω, τ ) −Z∞Zt−dτ0ds01 n −1s 2 exp(−s) ×Γ( n2 )pp× E d0 (x + 2 s(t − τ )Ω, τ )u(x + 2 s(t − τ )Ω, τ ) .

(2.2.74)88Несмещенной оценкой для Ku(x, t) будет случайная величинаs !n+ααα 2 Γ(n+α2 )g̃γη̃(x, t) = t 2x+2tϑΩ, x, t ×1αΓ( n2 )2s !n+α×u x+2 γtϑΩ, t − tϑ) +2r α 1n+ t 2 g̃2 x + 2 γtϑΩ, x, t − tϑ, t ×α2r n×u x+2 γtϑΩ, t − tϑ)) +2s !n+2+αα+12Γ()αn+2+α2+ t2h̃1 x + 2 γtϑΩ, x, t ×nαΓ( 2 )2s !n+2+αtϑΩ, t − tϑ) +×u x+2 γ2s !α nn+2+ t 2 × h̃2 x + 2 γtϑΩ, x, t − tϑ, tα2s !n+2tϑΩ, t − tϑ +×u x+2 γ2s !n+1)2Γ(1n+12+ t2d0 x + 2 γtδΩ, t − tδ A−1 (x, t)Ω ×nΓ( 2 )2s !n+1×u x+2 γtδΩ, t − tδ −2r n− td0 x + 2 γtθΩ, t − tθ ×2r n×u x+2 γtθΩ, t − tθ , (2.2.75)2где, как и формуле (2.2.53), случайные величины ϑ, δ, θ распределены на отрезке[0, 1], а γ(m) — случайная величина, имеющая гамма распределение с параметром m.