Решение систем линейных уравнений

Рассмотрим применение метода Монте-Карло к решению системы линейных уравнений:

. (9.13)

Пусть система (9.13) приведена к виду

(9.14)

где норма матрицы удовлетворяет условию . Тогда система (9.14) имеет единственное решение.

Приведем без доказательства схематическое описание метода Монте-Карло для решения системы линейных уравнений (9.14) (доказательство приведено в [7]).

Подберем такие множители v_ij, чтобы решения p_ij уравнений удовлетворяли условиям:

1)

2)

Положим

Будем трактовать число как вероятность перехода некоторой частицы из состояния S_i в состояние S_j. Всего состояний :

S₁, S₂, …, S_{n +1},

причем граничное состояние S_{n +1} = Г соответствует остановке частицы, так как , т.е. частица не может выйти из состояния S_{n +1}.

Матрица с элементами называется матрицей перехода состояний {S_i}:

(9.15)

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

(9.16)

Поставим в соответствие траектории T_i случайное число X_i:

, (9.17)

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

Теорема (доказательство в [7]). Математические ожидания случайных величин (9.17) равны корням системы (9.14): .

Сформулируем алгоритм решения системы линейных уравнений (9.13) методом Монте-Карло.

1. Привести систему (9.13) к виду (9.14):

(9.18)

2. Подобрать такие числа v_ij, чтобы решения p_ij уравнений удовлетворяли условиям:

1) 2) (9.19)

Вычислить -й столбец и -ю строку матрицы перехода:

(9.20)

3. Выполнить для каждого :

3.0. Присвоить значение x_i = 0;

3.1. Для каждого , (где N число испытаний случайной величины X_i), выполнить статистические испытания по алгоритму:

3.1.0. Присвоить m = 0; i_m = i; v = 1; x_i = x_i + β_i

3.1.1. Присвоим переменной очередное значение: m = m + 1;

3.1.2. Возьмем случайное число ξ из интервала (0; 1), и выберем значение i_m по следующему правилу:

3.1.3. Если j < n + 1, перейти к 3.1.1.

3.1.4. Для каждого вычислить

, ;

3.2. Вычислить x_i = x_i / N.

Приведем простой пример решения системы линейных уравнений методом Монте-Карло.

Пример 9.3. Решить систему уравнений

Решение.Приближенное решение x₁ = 0,68; x₂ = 0,93 получено методом итераций с точностью до 0,01.

Продемонстрируем применение алгоритма метода Монте-Карло, выполнив вычисления вручную. Запишем матрицу α и вектор β:

Очевидно, что можно выбрать числа p_ij и v_ij следующим образом:

Тогда выполняются условия (9.19). Вычислим 3-й столбец и 3-ю строку матрицы перехода по формулам (9.20):

Результаты дальнейших вычислений приведены в таблице 9.4, в которой показано, как вычисляется приближенное значение x₁. Мы используем здесь методику записи вычислений из учебника [7], прибегая к более подробному пояснению.

Проведем N = 10 испытаний случайной величины X₁. Для этого может понадобиться больше 10-ти случайных чисел. Вычислим с помощью программы Excel последовательность из 30-ти случайных чисел из интервала
(0; 1). Для этого введем в ячейку A1 формулу =СЛЧИС() и маркером заполнения протянем ячейку A1 до A30. Форматируем диапазон A1:A30 с помощью команды меню «Формат — Ячейки — Числовой — Число десятичных знаков 2». В первом столбце таблицы 9.4 записаны полученные в Excel случайные числа.

Так как мы вычисляем X₁, то начальное состояние S_i есть S₁. Всего состояний три: S₁, S₂ и S₃ = Г, Г — граничное состояние.

Согласно алгоритму, нам нужно определить в какой из следующих трех интервалов

[0; 0,4), [0,4; 0,4 + 0,1) = [0,4; 0,5), [0,5; 1)

попадает очередное случайное число ξ.

Если число попадает в первый интервал, то очередное состояние есть S₁, т.е. имеем переход S₁ → S₁. Если число попадает во второй интервал, то имеем переход S₁ → S₂, если в третий, то S₁ → S₃, т.е. частица переходит граничное состояние.

Если начальное состояние есть S₂, то надо проверять попадание частицы в интервалы

[0; 0,2), [0,2; 0,7), [0,7; 1).

Первое случайное число 0,07 попадает в интервал [0; 0,4), следовательно, частица переходит из S₁ в S₁. Второе случайное число 0,25 также попадает в интервал [0; 0,4), частица переходит из S₁ в S₁. Третье число 0,99 попадает в третий интервал [0,5; 1), поэтому частица переходит в граничное состояние S₃ = Г. Получаем траекторию S₁ → S₁ → S₁ →Г.

Теперь вычислим по формуле (9.17) значение X₁. Значения индексов i_m следующие: i₀ = 1, i₁ = 1, i₂ = 1, i₃ = 3 (m = 3).

.

Аналогично вычислим второе значение X₁. Следующее случайное число 0,73 попадает в третий интервал, получаем траекторию S₁ →Г и значение .

При вычислении третьего значения X₁ первое случайное число 0,14 попадает в первый интервал, второе случайное число 0,06 опять попадает в первый интервал, третье случайное число 0,54 попадает во второй интервал. Имеем часть траектории S₁ → S₁ → S₁ → S₂. Для состояния S₂ следующее случайное число 0,58 проверяем на интервалах

[0; 0,2), [0,2; 0,7), [0,7; 1).

Число 0,58 попадает во второй интервал, а следующее 0,68 — в третий. Траектория имеет вид S₁ → S₁ → S₁ → S₂→ S₂→Г и значение X₁ для этой траектории равно

Таблица 9.4

Аналогично можно вычислить среднее арифметическое для X₂.

Результаты вычислений приведены в таблице 9.5. Здесь начальное состояние есть S₂.

Таблица 9.5

Сравнивая полученные значения X₁ = 0,63; X₂ = 0,98 с приближенными корнями x₁ = 0,68; x₂ = 0,93, видим, что получены удовлетворительные результаты. Очевидно, для повышения точности необходимо провести значительно большее число испытаний.