Ім'я файлу: 1.pdf
Розширення: pdf
Розмір: 1064кб.
Дата: 03.04.2020
скачати

f (x) =
2 3
(x + 1). Запишем интеграл так
I =
1
Z
0
e x
2 3
(x + 1)
·
2 3
(x + 1)dx.
Таким образом, интеграл I представлен в виде математического ожидания функции e
x
2 3
(x+1)
. В качестве искомой оценки примем выборочную среднюю
(для простоты ограничимся десятью испытаниями):
1 10 10
X
i=1
e x
i
2 3
(x i
+ 1)
= 0.15 10
X
i=1
e x
i x
i
+ 1
,
где x i
– возможные значения X, которые надо разыграть по известной плот- ности f (x) =
2 3
· (x + 1).
Имеем
2 3
x i
R
0
(x + 1)dx = r i
. Отсюда находим явную формулу для разыгрыва- ния возможных значений X:
x i
=
√
1 + 3 · r i
− 1.
Результаты 10 испытаний приведены в следующей таблице:
Номер испытания, i
1 2
3 4
5 6
7 8
9 10
r i
0.100 0.973 0.253 0.376 0.520 0.135 0.863 0.467 0.354 0.876
x i
=
√
1 + 3r i
− 1 0.140 0.980 0.326 0.459 0.600 0.185 0.894 0.550 0.436 0.905
e x
i
1.150 2.664 1.385 1.582 1.822 1.203 2.445 1.773 1.546 2.472
x i
+ 1 1.140 1.980 1.326 1.459 1.600 1.185 1.982 1.550 1.436 1.905
e xi x
i
+1 1.009 1.345 1.044 1.084 1.139 1.015 1.291 1.118 1.077 1.298
Сложив числа последней строки таблицы, получим
10
P
i=1
e xi x
i
+1
= 11.42. Та- ким образом, искомая оценка равна I
∗
2
= 0.15 · 11.42 = 1.713. Заметим, что точное значение интеграла I = 1, 718...
27

1.4.3
Способ, основанный на истолковании интеграла как пло- щади
Вычислим интеграл (1.28), где подынтегральная функция неотрицательна и ограничена
0 ≤ ϕ(x) ≤ c,
исходя из истолкования интеграла как площади. Введем в рассмотрение двумерную случайную величину (X, Y ), распределенную равномерно в прямоугольнике D с основанием (b−a) и высотой c , плотность вероятности которой f (x, y) =

1
(b−a)c для точек, принадлежащих D;
0,
вне D.
Составляющая X распределена в интервале (a, b) равномерно с плотно- стью распределения
1
b−a
; составляющая Y распределена в интервале (0, c)
с плотностью
1
c
. Если разыграно N точек (x i
, y i
), принадлежащих прямо- угольнику D , из которых n точек оказались под кривой, то отношение площади, определяемой интегралом (1.28), к площади прямоугольника D
b
R
a
ϕ(x)dx
(b − a)c
≈
n
N
(1.39)
Отсюда b
Z
a
ϕ(x)dx ≈ (b − a)c n
N
(1.40)
Таким образом, в качестве оценки интеграла (1.28) можно взять
I
∗
3
≈ (b − a)c n
N
,
(1.41)
где N – общее число случайных точек (x i
, y i
), принадлежащих D ; n –
число случайных точек, которые расположены под кривой y = ϕ(x).
Пример 1.4
Найти оценку интеграла
2
R
0
(4 − x
2
)dx.
28

Решение.
В интервале (0, 2) подынтегральная функция ϕ(x) = 4 − x
2
неотрицательна и ограничена, причем ϕ(x) ≥ ϕ(0) = 4; следовательно,
можно принять c = 4.
Введем в рассмотрение двумерную случайную величину (X, Y ), распре- деленную равномерно в прямоугольнике D с основанием b − a = 2 и вы- сотой c = 4, плотность вероятности которой f (x, y) =
1 2·4
=
1 8
. Разыграем
N = 10 случайных точек (x i
, y i
), принадлежащих прямоугольнику D. Учи- тывая, что составляющая X в интервале (0, 2) распределена равномерно с плотностью f
1
(x) =
1 2
и составляющая Y в интервале (0, 4) распределе- на равномерно с плотностью f
2
(y) =
1 4
, разыграем координаты случайной точки (x i
, y i
), принадлежащей прямоугольнику D, по паре независимых случайных чисел (r i
, R
i
):
1 2
·
x i
Z
0
dx = r i
,
1 4
·
y i
Z
0
dy = R
i
Отсюда x i
= 2 · r i
, y i
= 4 · R
i
. Если окажется, что y i
< 4 − x
2
i
, то точка
(x i
, y i
) лежит под кривой Y = ϕ(X) и в счетчик n надо добавить единицу.
Результаты десяти испытаний приведены в таблице.
Номер испытания, i r
i x
i
= 2 · r i
x
2
i
4 − x
2
i
R
i y
i
= 4 · R
i y
i
< 4 − x
2
i
1 0.100 0.200 0.040 3.960 0.973 3.892 1
2 0.253 0.506 0.256 3.744 0.376 1.504 1
3 0.520 1.040 1.082 2.918 0.135 0.540 1
4 0.863 1.726 2.979 1.021 0.467 1.868 5
0.354 0.708 0.501 3.499 0.876 3.504 6
0.809 1.618 2.618 1.382 0.590 2.360 7
0.911 1.811 3.320 0.680 0.737 2.948 8
0.542 1.084 1.175 2.825 0.048 0.192 1
9 0.056 0.112 0.013 3.987 0.489 1.956 1
10 0.474 0.948 0.899 3.101 0.296 1.184 1
Из таблицы находим n = 6. Искомая оценка интеграла I
∗
3
= (b−a)·c·
n
N
=
2 · 4 ·
6 10
= 4.8.
29

1.4.4
Способ выделения главной части
Пусть требуется вычислить интеграл (1.28). Введем в рассмотрение слу- чайную величину X, распределенную равномерно в интервале интегриро- вания (a, b) с плотностью f (x) =
1
b−a
. Допустим, что найдена такая функ- ция g(x), которая мало отличается от ϕ(x) и интеграл от которой мож- но вычислить, не прибегая к методу Монте-Карло. Тогда математическое ожидание функции
F (X) = (b − a) [ϕ(X) − g(X)] +
b
Z
a g(x)dx
(1.42)
равно искомому интегралу I.
Действительно, учитывая, что величина X распределена в интервале интегрирования (a, b) равномерно с плотностью
1
b−a
, получим
M [F (X)] = (b − a)
b
Z
a
[ϕ(X) − g(X)]
1
b − a dx +
b
Z
a g(x)dx =
b
Z
a
ϕ(x)dx = I.
Таким образом, в качестве оценки математического ожидания M [F (X)], а следовательно и интеграла I, можно принять среднее арифметическое N
значений функции F (X) :
I
∗
4
=
b − a
N
N
X
i=1
[ϕ(x i
) − g(x i
)] +
b
Z
a g(x)dx,
(1.43)
где x i
– возможные значения X, которые разыгрывают по формуле x
i
= (b − a) · r i
Пример 1.5
Найти оценку I
∗
4
интеграла
1
R
0
√
1 + x
2
dx.
Решение.
Так как
√
1 + x
2
= 1 +
x
2 2
+ ...
(|x| ≤ 1), то примем g(x) = 1 +
x
2 2
. Тогда, полагая число испытаний N = 10, имеем оценку
I
∗
4
=
1 10 10
X
i=1
q
1 + x
2
i
−

1 +
x
2
i
2

+
1
Z
0

1 +
x
2 2

dx.
30

Выполнив элементарные преобразования, получим
I
∗
4
=
1 20 10
X
i=1

2
q
1 + x
2
i
− x
2
i

+
1 6
Учитывая, что a = 0, b = 1, возможные значения x разыграем по формуле x
i
= a + (b − a) · r i
= r i
Результаты вычислений приведены в таблице.
Номер испытания, i
1 2
3 4
5 6
7 8
9 10
x i
= r i
0.100 0.973 0.253 0.376 0.520 0.135 0.863 0.467 0.354 0.876
x
2
i
0.010 0.947 0.064 0.141 0.270 0.018 0.745 0.218 0.125 0.767 1 + x
2
i
1.010 1.947 1.064 1.141 1.270 1.018 1.745 1.218 1.125 1.767
p1 + x
2
i
1.005 1.395 1.032 1.068 1.127 1.009 1.321 1.104 1.061 1.329 2
p1 + x
2
i
− x
2
i
2.000 1.843 2.000 1.995 1.984 2.000 1.897 1.990 1.997 1.891
Сложив числа последней строки таблицы, найдем сумму 19.597. Таким образом, получим искомую оценку интеграла I
∗
4
=
19,597 20
+
1 6
= 1.145. Заме- тим, что точное значение I = 1.147.
31

1.4.5
Вычисление кратных интегралов
Пусть функция y = f (x
1
, ..., x m
) непрерывна в ограниченной замкнутой области G и требуется вычислить m-кратный интеграл I по области G:
I =
Z
· · ·
Z
G
f (x
1
, . . . , x m
)dx
1
. . . dx m
(1.44)
Геометрически число I представляет собой (m + 1)-мерный объем верти- кального цилиндрического тела в пространстве Ox
1
x
2
. . . x m
y, построенно- го на основании G и ограниченного сверху данной поверхностью y = f (x),
где x = (x
1
, . . . , x m
).
Преобразуем интеграл так, чтобы новая область интегрирования Ω це- ликом содержалась внутри единичного m-мерного куба. Пусть область ин- тегрирования G расположена в m-мерном параллелепипеде a k
≤ x k
≤ b k
(k = 1, . . . , m).
Сделаем замену переменных x
k
= a k
+ (b k
− a k
)ξ
k
(k = 1, . . . , m).
(1.45)
Тогда m-мерный параллелепипед преобразуется в m-мерный единичный куб 0 ≤ ξ
k
≤ 1(k = 1, . . . , m), и, следовательно, новая область интегриро- вания Ω будет целиком расположена внутри этого единичного куба.
Вычислим якобиан преобразования:
D(x
1
, . . . , x m
)
D(ξ
1
, . . . , ξ
m
)
=
b
1
− a
1 0
0 0
b
2
− a
2 0
0 0
. . . b m
− a m
= (b
1
−a
1
)·(b
2
−a
2
)·. . .·(b m
−a m
).
(1.46)
Таким образом,
I = (b
1
− a
1
) · (b
2
− a
2
) · . . . · (b m
− a m
) · J,
(1.47)
где
J =
Z
· · ·
Z
Ω
f [a
1
+ (b
1
− a
1
)ξ
1
, . . . , a m
+ (b m
− a m
)ξ
m
]dξ
1
. . . dξ
m
(1.48)
и область интегрирования Ω содержится внутри m-мерного единичного ку- ба
0 ≤ ξ
k
≤ 1
(k = 1, . . . , m).
(1.49)
32

Пусть
F (ξ
1
, . . . , ξ
m
) = f [a
1
+ (b
1
− a
1
)ξ
1
, . . . , a m
+ (b m
− a m
)ξ
m
].
(1.50)
Выберем N равномерно распределенных на отрезке [0, 1] последователь- ностей случайных чисел:
ξ
(1)
1
,
ξ
(1)
2
, . . . , ξ
(1)
m
;
ξ
(2)
1
,
ξ
(2)
2
, . . . , ξ
(2)
m
;
ξ
(N )
1
, ξ
(N )
2
, . . . , ξ
(N )
m
(1.51)
Рассмотрим случайные точки M
i
(ξ
(i)
1
, ξ
(i)
2
, . . . , ξ
(i)
m
). Выбрав достаточно большое N число точек M
1
, . . . , M
N
проверяем, какие из них принадлежат области Ω. Пусть n точек принадлежат области Ω.
Заметим, что относительно границы области Ω следует заранее догово- риться, причисляются ли граничные точки, или часть их, к области Ω, или не причисляются к ней. В общем случае при гладкой границе это не име- ет существенного значения; в отдельных случаях нужно решать вопрос с учетом конкретной обстановки.
Взяв достаточно большое n точек M
i
, F
ср приближенно можно положить равным:
e
F
ср
=
1
n n
X
i=1
F (M
i
);
(1.52)
откуда оценка интеграла J выражается формулой e
J = ˜
F
ср
· Ω,
(1.53)
где под Ω понимается m-мерный объем области интегрирования Ω. Если вычисление объема Ω затруднительно, то можно принять e
Ω =
n
N
(из опре- деления геометрической вероятности).
Таким образом, имеем оценку e
I искомого интеграла I:
e
I = V · e
F
ср
· e
Ω,
(1.54)
где V = (b
1
−a
1
)·. . .·(b m
−a m
) – объем параллелепипеда, ограничивающего исходную область интегрирования G. Отсюда получаем e
I =
V
N
n
X
i=1
F (M
i
).
(1.55)
33

В частном случае, когда Ω есть единичный куб, проверка принадлежно- сти точек M
i области Ω становится излишней, то есть n = N и мы имеем e
I =
V
N
N
X
i=1
F (M
i
).
(1.56)
Пример 1.6
Вычислить интеграл
1
R
0
dx
1
R
x
(x + y)dy. Произвести 10 испытаний.
Решение.
Область интегрирования ограничена линиями y = x, x = 1,
x = 0 и, очевидно, принадлежит единичному квадрату. Площадь области интегрирования (прямоугольного треугольника)
S =
1 · 1 2
= 0.5.
Используем формулу
I
∗
1
= S ·
n
P
i=1
f (x i
, y i
)
n
,
где n – число случайных точек (x i
, y i
), которые принадлежат области интегрирования; у этих точек y i
> x i
(при каждом испытании, в котором это условие выполняется в счетчик n записывают единицу). Пары незави- симых случайных чисел (x i
, y i
) берем из таблицы приложения 1, начиная с первой строки сверху. Результаты 10 испытаний приведены в следующей таблице.
34

Номер испытания, i x
i y
i
Счетчик y i
≥ x i
f (x i
, y i
) = x i
+ y i
1 0.100 0.973 0.253 0.376 2
0.100 0.973 1
1.073 3
0.520 0.135 4
0.863 0.467 5
0.354 0.876 1
1.230 6
0.809 0.590 7
0.911 0.737 8
0.542 0.048 9
0.056 0.489 1
0.545 10 0.474 0.296
P
4 3.477
Из таблицы находим n = 4,
4
P
i=1
f (x i
, y i
) = 3.477. Подставив эти числа в формулу, получим искомую оценку:
I
∗
1
= 0.5 ·
3.477 4
= 0.435.
Сравнительно большое расхождение полученной оценки с точным зна- чением I = 0.5 объясняется малым числом испытаний.
35

1.5
Погрешность интегрирования методом Монте-Карло
1.5.1
Формула оценки погрешности интегрирования методом
Монте-Карло
Итак, для оценки интеграла
I =
Z
· · ·
Z
G
f (x
1
, . . . , x m
)dx
1
dx
2
. . . dx m
имеем формулу e
I = V · e
F
ср
· e
Ω,
где V – объем параллелепипеда, ограничивающего область интегрирова- ния G;
e
Ω – объем новой области интегрирования Ω, полученной после замены пе- ременных (1.45);
e
F
ср
– оценка среднего значения подынтегральной функции, полученной по- сле замены переменных (1.45). Оценим погрешность интегрирования.
Из теории погрешностей известна формула оценки погрешности произ- ведения двух величин:
∆I = ∆(V · e
F
ср
· e
Ω) = V ·

∆F
ср
· e
Ω + ∆Ω ·
e
F
ср

(1.57)
По центральной предельной теореме
P
(
F
ср
− e
F
ср
≤ t
β
r
σ
2 1
n
)
≈ β;
P
(
Ω − e
Ω
≤ t
β
r
σ
2 2
N
)
≈ β;
т.е. с вероятностью β
∆F
ср
= t
β
r
σ
2 1
n
, ∆Ω = t
β
r
σ
2 2
N
;
при этом
36

e
F
ср
=
1
n n
X
i=1
F (M
i
), e
Ω =
n
N
Подставим значения погрешностей в формулу (1.57).
∆I = V ·
t
β
r
σ
2 1
n
· e
Ω + t
β
r
σ
2 2
N
·
e
F
ср
!
(1.58)
Выполнив элементарные преобразования, при известных значениях σ
1
и
σ
2
, для оценки погрешности имеем формулу
∆I = V · t
β


σ
1
·
e
Ω
√
n
+ σ
2
·
e
F
ср
√
N


(1.59)
Однако, обычно при решении конкретных задач величины σ
2 1
и σ
2 2
ока- зываются неизвестными. Тогда их приближенные значения вычисляют по формулам
S
2 1
= ( e
F
2
)
ср
−
 e
F
ср

2
,
S
2 2
= e
Ω · (1 − e
Ω),
где e
F
ср
=
1
n n
X
i=1
F (M
i
), ( e
F
2
)
ср
=
1
n n
X
i=1
F
2
(M
i
), e
Ω =
n
N
Тогда для оценки погрешности имеем формулу
∆I = V · t
β


S
1
·
e
Ω
n
+ S
2
·
e
F
ср
√
N


(1.60)
Так как из того, что e
Ω =
n
N
следует, что n = e
Ω · N, то
∆I = V · t
β


S
1
p e
Ω
√
N
+
S
2
e
F
ср
√
N


(1.61)
При достаточно большом N , а следовательно и n, e
F
ср
→ F
ср и e
Ω → Ω (а
F
ср и Ω не зависят от N ). Предположим, что величина
V · t
β

S
1
p e
Ω +
e
F
ср
S
2

37

с ростом N изменяется незначительно (принимает приблизительно одно и то же значение). Обозначим c = V · t
β

S
1
p e
Ω +
e
F
ср
S
2

Тогда погрешность ∆ =
c
√
N
. Если полученная погрешность ∆ велика, то мы назначаем желаемую погрешность (точность) ∆
0
.Тогда количество ис- пытаний N
0
, необходимое для достижения указанной точности приближен- но можно положить равным N
0
=
c
∆
0


2
и вычислить значение интеграла с новым количеством испытаний.
38

1.5.2
Некоторые способы уменьшения дисперсии
Количество испытаний, необходимое для вычисления интеграла с задан- ной точностью, зависит, как было показано выше, от дисперсии соответ- ствующей случайной величины. Существуют различные приемы преобра- зования задачи, позволяющие уменьшать дисперсию и, следовательно, вре- мя вычислений.
Если часть задачи можно решить аналитически, то, используя это частичное решение, обычно удается построить метод Монте-Карло для решения всей задачи с меньшей дисперсией. Правда, вообще говоря,
построенные таким путем методы могут оказаться более трудоемкими и в конечном счете невыгодными. Мы будем говорить об аналитическом интегрировании, хотя в некоторых случаях такую же роль может сыграть численное интегрирование, если точность этого интегрирования значи- тельно выше, чем точность метода Монте-Карло.
Выделение главной части
Если главную часть задачи можно вычислить аналитически, то, как правило, выгодно считать методом Монте-Карло не всю задачу, а только
«поправку» — разницу между всей задачей и главной частью. Уменьшение дисперсии при этом может оказаться очень значительным.
Интегрирование по части области
Допустим, что мы умеем (аналитически) вычислить интегралы по некоторой части области G. Тогда выгодно представить интеграл в виде суммы интегралов по областям B и G − B, и методом Монте-Карло вычислить только интеграл по области G − B. Чем область B больше, тем заметнее будет понижение дисперсии.
Интегрирование по части переменных (понижение порядка интеграла)
Если аналитически взять интеграл по некоторым из переменных, а по остальным переменным использовать тот же метод Монте-Карло, то дис- персия уменьшится. Правда, нередко бывает, что после интегрирования по некоторым из переменных получаются более сложные формулы счета и,
несмотря на уменьшение дисперсии, трудоемкость возрастает.
39

Глава 2
Численный эксперимент
2.1
Алгоритм вычисления интегралов методом Монте-Карло

1   2   3   4   5

скачати