Соболь И.М. Численные методы Монте-Карло (1973) (1186217), страница 32
Текст из файла (страница 32)
значения величины ьЩ не изменят. Если матрица А содержит целую строку из нулей: а„,=ааа=„.. =аа„=О, то, конечно, нельзя все соответствующие вероятности выбрать нулевыми: это противоречило бы условию (60). Тогда можно положить Ран= 1, !Таа 0 ПРИ РЧЬСЕ, так ЧТО 7!аз=бар. 5.3. Вычисление одной компоненты решения. Иногда встречаются задачи, сводящиеся к системе (55), где, однако, нас интересует не все решение г= (гь ..., г ), а только одна из неизвестных, например г„. Методы Монте-Карло позволяют приближенно оценить одну эту компоненту. Для этого достаточно в качестве чр выбрать единичный вектор е'м= (О,..., О, 1, О,..., 0), в котором лишь на и-м месте стоит 1. Тогда скалярное произведение равно (тр,г) = (е<'>,г) = ~ ен7г = г,.
а В качестве начальных вероятностей можно выбрать ра= е„" или, другими словами, начинать цепь с Йа=г. (г! !гл а РЕШЕНИЕ ЛННЕННЫХ УРЛВНЕННЙ 196 Итак, строим цепи г-«й!-»)га-»... -«Й!-»..., вы- числяем вдоль цепей веса нь ма," а!,а, о л ... а, (66) игарка,а» " ра!,а, и случайные величины ь (е"! = Х Ото»п ,=е Если количество цепей )у' достаточно велико, то г, — ~ч Ь (е<г!)„ зси (67) где ь[е!"!~),— значение ь[е!'!), полученной на з-й цепи.
БА. Обрашение матрицы. Пусть задана квадратная матрица В= ((! р), 1~се, Р(гн. Требуется вычислить элементы обратной матрицы С=(с а), удовлетворяющей соотношению ВС=Е, где Š— единичная матрица: О 1 ... О Если все собственные значения р„ матрицы В удовлетворяют условию (р.— 1! <1, то задача эта может быть решена методом последовательных приближений* ). В самом деле, обозначим через А матрицу А =Š— В. (68) Легко проверить, что собственные значения матрицы А равны ! — р„и, следовательно, по абсолютной величине меньше единицы.
В силу этого сходится ряд С=~~Р, А! и ~=О ОЭ 0 ВС= (Š— А) С= С АС ~2~ А,' — 2~ РА =Е **). =о *) Более оби!ий случай см. упражнение 6 стр. 209. ") Для запоминания можно заметить, что рял (Š— А) ОР « ° = ~Ч~~А' аиало ичен геометричесаой прогрессии (1 — 4) '= ~~Р~ о' ° г-о !-е $ и Решение линепных АлгеБРлическнх систем 199 так что его г-я компонента есть В частном случае, когда [= ео ~, получим, что г„= с,тч что и требовалось доказать. Подставив в (бу) вместо 1 вектор ео 1, получим со. ответствующий этому случаю метод Монте-Карло: при достаточно большом Л' с„' ж —,, ~~э~ ~ ()те (69) здесь снова з — номер цепи, а сумма с условием (1)м,=г') означает, что суммируются толькО те ((Уи для которых номер Й,=г'.
бм.1. Численный пример: нанти матрицу С, обратную к матрице В= — 2 1!2 О собственные аначеиии которой, очевидно, равны р, ра 1/2, ра !. Сперва по формуле (ба) вычислим матрицу А= 2 1/20 А(окажем, что если матрица А удовлетворяет условиям теоремы 8, то для того, чтобы вычислить элемент с„,' матрицы С, можно использовать цепи п. 5.3; надо лишь, кроме тр=е"', выбрать еще [=ео'1. В самом деле, мы видели, что М~[ео11=ам где г,— г-я компонента решения уравнения (бб). Однако это решение равно РЕШЕ))НЕ ЛИНЕЙНЫХ УРАВНЕНИЙ !Гл 5 200 Затеи выбереч матрицу вероятностей перехода, допустимую по от ношению к матрице А, например: 1 0 0 1 1 (р„р) = —,, —,, о 1 0 О (70) Веса вдоль такой пепи меняются на множитель оп)рп 1/з 1, 1 йг =1, (Р = —, 1Г.= —,...
2' ' 4' "' Осредиение в (69) не нужно, достаточно одной цепи для получении точного результата 1 1 сз)= г )з')=1+ 2+ 4+... 2, зюя Цз)=Н сзз — — ~ йг) — — 0; !з !51=а) 2) Столь же просто оцениваются элементы сзь сзз н сзз' цепь начинается с йа г=31 так как в третьей строке матрицы (70) лишь ОДИП НЕНУЛЕВОЙ ЭЛЕМЕит Рм = 1, тО йз 1, ДаЛЕЕ Аз=1, йз= 1, ... Итак, все цепи будут одинаковыми йа=З, й 1, й =1, йз=(,..., а соответствующие веса (с учетом того, что аз)/раз=1) равны 1, 1 йуа=! )" з=! йгз= 2' з 4 Поэте))у 1 ! сзз= ~и~~~ ~йтг=! ) 2+ 4+...=2з 11Щ)=1) с = Х )р)=0, !Оа,=г) 3) И толыго для вычисления см, сзз и саз нужно действительно строить случайные цепи, начинающиеся с да=)=2, и производить расчеты.
Впрочем, в этом случае различные типы цепей легко классифицируются, ибо во второй строке матрицы (70) лишь рз) и рзз отличны от нуля. 1) Чтобы оценить элементы сзь г„ и с,з, надо строить цена, начинающиеся с йа=г 1 Таи как в первой строке л)атрицы (70) лишь рпФО, тс есе цепи окажутся одинаковыми: да=), 4 з1 Решение линеиных АлГеБРАических систез! 291 Цепь, начала которой до=2, может несколько раз оставаться в й~ — — йз=...=йг — — 2. Ио если окажюся, что й, ! — — 1, то дальше !+ уже все йг+! = а~+э=...=!.
Следовательно, общий впд такой цепи (О<!<оо) йо = йт = .. = Аг = 2, йо+! = Фю э — — ... —— 1. (т!) Соответствующие такой цепи веса легко вычпслиток так как ам(ры =1, ам)рм=4, апlрп='!о, то )ио = )Угз = " — — (Ггг — — 11 1о'!+! = 4; В'г+э — — 2, йн рз = 1... ° Иитересующпе иас суммы весов вдоль ~акой цепи равны йг,,=4+2+!+ 1)2+... =8, яэ,-!) Иг=т+ 1, ~чг, !Р,.=О. (г(э!=э) (Иэг=з) Из этих соотношений сразу следует, что си=8, сто=о, Оставшийся элемеиг с,о можио вычислить аналитически, если записать, что вероятность получигь цепь (7!) равна (см: вача.ю стр. 39) (Ро„)' Ро, = (1г2)Г+~, и вместо формулы (69) воспользоваться точным выражением длп математического ожидания с о = м ~чР~ и' = ~и~ (! -)-1) 2 <г+з! =2. (где!А=» г-о Итак, обратная матрица С равна С= 8 2 О Конечно, в общем случае (когда в матрице А много ненулевых элемеитов) разнообразие цепей столь велико, что заменить Грормулу (89) аналитическим расчетом невозможно.
Заметим также, что по цепям, иачииающимся с йо=г, можно одновременно вычислять все сгг' с г'= 1, 2, ..., т, ибо каждое Ф'г входит в сумму яэ! г') при одном и только одном г'. 5.5. Решение дифференциальных уравнений Лапласа и Пуассона. В ограниченной связной области сг плоскости х, у с простой границей бо рассмотрим РЕШЕНИЕ ЛИНЕИНЫХ УРАВНЕНИИ 1гл. 5 дифференциальное уравнение с частными производными дхи/дхз+д'и/дуз=с (х, у), (72) где и=и(х, у) — искомая функция.
Уравнение (72) при г (х, у) =О называется уравнением Лапласа, а при г (х, у) чйΠ— уравнением Пуассона. Предположим, что на границе 6и задана некоторая функция у(х, у) (часто пишут К(5), где Я вЂ” длина дуги границы, отсчитываемая от какой-нибудь фиксированной точки). Требуется найти такое решение и(х, у) уравнения (72), которое на границе совпадает с у(х, у): и(х, у) в,=д. (73) Задачу об отыскании решения уравнения (72), удовлетворяющего грал ничному условию (73), называют задачей Дирихле.
Для приближенного решения этой задачи [63, 88] выбирают на плоскости достаточно мелкую -гиеииииии глм 1-Г'и-Г которое можно переписать в инде и) ~.— (1/4) (и ~ ~+ и „, + и. ~ + и.1+ — )РР ~). (74) квадратную сетку с шаРис. 53. гом й (рис. 53). Координаты узлов этой сетки пусть будут х,=/11, у~=1Й, а значения и(х„у,) н г(х„у) для краткости обозначим и,, и г, . Узел (/, 1) называют внутренним, если и он, и все четыре соседних с пим узла (1 — 1, 1), (1+1, 1), (1, 1 — 1), (1, 1+!) принадлежат 6+6', в противном случае узел (1, 1), принадлежащей 6+6з, называют граничныи. Во внутреннем узле (хь у,) уравнение (72) заменяется разностным уравнением и+ 1 — 2и)~+ и), и 1+1 — 2и; + и ь + ни 2 2! Решение линейных ллгеБРлических систем 2ОЗ В граничных узлах полагают (75) Нз ~=й 1. (Значения д1,! «Сносят» с ближайших точек границы 6'.) Решение алгебраической системы (74) — (75) прн Й- О приближается к решению задачи Дирихле для уравнения (?2).
Если перенумеровать все узлы, принадлежащие О+!2' (в произвольном порядке), и переписать в том же порядке уравнения (74) и (75), то получив! систему вида (55) на= Х пазнз+/а, с!=1 2, °,, 1н а=! '(п2 — количество узлов), с весьма своеобразной матрицей А: внутреннему узлу с номером а отвечает строка аа1,..., а,, в которой 4 элемента равны 1/4, а остальные — нули; граничному узлу с номером с! отвечает строка а 1 — — а 2= ... =а„„=О; все диагональные элементы а а=О.
(Можно доказать, что все собственные значения такой матрицы по абсолютной величине меньше сд1шицы.) Свободные члены этой системы равны = — (1/4)ЙЕРа, если узел номер а внутренний, и / =Д, если узел номер а граничный. Воспользуемся методом п. 5.3. и построим метод Монте-Карло для расчета и, — значения решения в одном (заранее заданном) узле. Выберем матрицу переходов паз, если узел номер а внутренний, Раз = баз, если узел помер с! граничный ',(б з — снмвол Кронекера: б =1, б„з=О при аФй).
Процесс построения цепи по такому закону оказывается очень наглядным: 1) начинаем с йч —— г; 2) если узел й, внутренний, то с одинаковой вероятностью '/, выбираем в качестве й1+! номер одного из соседних с ним узлов; 3) если узел й! граничный, то цепь останавливается: й =/2!+1 = й!+2 — ° ° ° Расчет весов 1Р', вдоль такой цепи также чрезвычайно прост: пока цепь не попала на границу )Р',= В'1=... „.=%'1=1; далее (Р', 1= 1Р, 2=...
=О. Рншпнив лшчегпчых уи»внвнин (гл я Позтому случайная величина (67) оказывается равной ' = /», + /», +... + /»г (76) где й; — номер первого выхода цепи на границу. В (76) все /»„..., /, вычисляются по формуле », — (1/4)йтг„и лишь последнее /», равно значению ам 3 а и е ч а н и е. Если вместо граничных условий (73)' заданы более сложные условия, например: *(с1и+ся(ди/дх)+са (ди/ду)1о =д„ тг» уравнения (76) наряду с иь, будут содержать также значения иа в некоторых соседних узлах. И случайная цепь, попав на границу, останавл: виться не будет.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
