Приближенное вычисление интегралов
Пусть требуется вычислить определенный интеграл
На этом занятии мы познакомимся с простыми способами приближенного вычисления: формулой прямоугольников, формулой трапеций, формулой Симпсона или параболическим интегрированием, методом Монте-Карло.
3.1. Формула прямоугольников
Будем исходить из геометрических соображений и рассматривать определенный интеграл
Рис. 48
Идея, которая привела к понятию определенного интеграла заключалась в следующем. Разбить всю фигуру на полоски одинаковой ширины dx = (b - a)/n, а затем каждую полоску заменить прямоугольником, высота которого равно какой-либо ординате (см. рис. 48).
Тогда получится следующая формула:
где xi
<= ci <= xi+1 (i = 0, 1, ..., n-1). Площадь криволинейной фигуры заменится площадью сумм прямоугольников. Эта приближенная формула и называется формулой прямоугольников.
Практически, в качестве точки ci берут середину промежутка [xi , xi+1], т. е.
Нетрудно составить процедуру вычисления такой суммы.
{ Вычисление интеграла методом прямоугольников }
{ Rectangle - прямоугольник }
Procedure Rectangle(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a + dx/2;
f := fx(c);
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
Теперь возникает вопрос о числе точек деления - n, который напрямую связан с точностью вычисления интеграла.
Известно, что точность вычисления интеграла по формуле прямоугольников оценивается по дополнительному члену, который выражается формулой:
значит оценивать точность вычисления можно по модулю этого остаточного или дополнительного члена, |Rn|.
Вторую производную вычислить мы сможем с достаточно высокой степенью точности, применив следующую функцию:
{ Функция вычисления второй производной }
Function
derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
Осталось выяснить, в какой точке промежутка интегрирования [a, b] находить значение этой производной. Остаточный член не требует строго определенного значения аргумента из этого промежутка, поэтому можно выбрать любое значение, не выходящее за пределы интервала [a, b]. Сразу возникает мысль вычислить вторую производную в середине промежутка, т.е. в точке (a + b)/2. Но представьте себе ситуацию, когда промежуток [-1, 1] или [-6.28; 6.28], тогда середина этого отрезка - точка 0 и значение производной будет равно нулю, а значит для числа точек деления n может быть установлено значение любое, даже 1, что, конечно, не даст требуемой точности вычисления интеграла.
Итак, следующая проблема, в какой точке промежутка находить значение производной?
Можно найти наибольшее значение производной на промежутке интегрирования [a, b]. Это можно сделать с помощью процедуры:
{ Определение наибольшего значения второй производной }
Procedure Maximum(a, b, eps : real; var
max : real);
var
dx, x : real;
begin
dx := 0.1; x := a;
max := abs(derivat2(x, eps));
while x<= b do
begin
x := x + dx;
if max < abs(derivat2(x, eps))
then max := abs(derivat2(x, eps))
end
end;
В ней установлен шаг - 0.1 и затем, в каждой следующей точке промежутка определяется значение второй производной и находится наибольшее.
Теперь составим процедуру определения числа точек деления, а затем и полностью программу.
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps, max : real; var
n : integer);
var
d : real;
begin
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while (max*d)/(24*n*n) >= eps do n := n+1;
end;
Программа вычисления интеграла по формуле прямоугольников
{ Вычисление интеграла по формуле прямоугольников }
Program Jntegral_Rectangle2;
uses WinCrt;
var
a, b, eps, j : real;
n : integer;
{----------------------------------------------------------------------------------------}
Function fx(x : real) : real;
begin
fx := sin(x)
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления второй производной }
Function derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
{----------------------------------------------------------------------------------------}
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps : real; var
n : integer);
var
dy2, d, c, dx : real;
begin
c := (a + b)/2;
dy2 := derivat2(c, eps);
if dy2 = 0
then
begin
c := a; dx := (b - a)/10;
while derivat2(c, eps) = 0 do c := c + dx;
dy2 := derivat2(c, eps)
end;
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while abs(dy2*d)/(24*n*n) >= eps do n := n+1;
end;
{----------------------------------------------------------------------------------------}
{ Вычисление интеграла методом прямоугольников }
Procedure Rectangle(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a + dx/2;
f := fx(c);
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления интеграла '); readln(eps);
Number(a, b, eps, n); Rectangle(a, b, n, j);
writeln('Значение интеграла равно ', j:6:t(eps));
writeln('С точностью до ', eps:1:t(eps))
end.
Такой способ можно использовать, но он может быть долгим, особенно, если большой промежуток интегрирования, малый шаг и более высокая точность вычисления производной.
Для определения числа точек деления можно применить и другой прием.
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps : real; var
n : integer);
var
dy2, d, c, dx : real;
begin
c := (a + b)/2;
dy2 := derivat2(c, eps);
dx := (b - a)/10;
if dy2 = 0
then
begin
c := a;
while derivat2(c, eps) = 0 do c := c + dx;
dy2 := derivat2(c, eps)
end;
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while abs(dy2*d)/(24*n*n) >= eps do n := n+1;
end;
Как работает эта процедура? Вычисляется значение второй производной в середине промежутка: c := (a + b)/2. Если она не равна нулю, тогда все в порядке, вычисляется в зависимости от заданной точности число точек деления n. Если значение второй производной равно нулю в середине промежутка интегрирования: dy2 = 0, тогда устанавливается шаг dx := (b - a)/10 и начинается цикл, в котором вычисляется производная через каждый промежуток dx, начиная от точки a, пока производная равна нулю цикл продолжается и заканчивается как только она не станет равной нулю.
Такой прием экономит время работы программы, производную можно вычислить с большей точностью, да и сама программа становится менее громоздкой:
Program Jntegral_Rectangle2;
uses WinCrt;
var
a, b, eps, j : real;
n : integer;
Function fx(x : real) : real;
begin
fx := sin(x)
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления второй производной }
Function derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
{----------------------------------------------------------------------------------------}
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps : real; var
n : integer);
var
dy2, d, c, dx : real;
begin
c := (a + b)/2;
dy2 := derivat2(c, eps);
if dy2 = 0
then
begin
c := a; dx := (b - a)/10;
while derivat2(c, eps) = 0 do c := c + dx;
dy2 := derivat2(c, eps)
end;
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while abs(dy2*d)/(24*n*n) >= eps do n := n+1;
end;
{----------------------------------------------------------------------------------------}
{ Вычисление интеграла методом прямоугольников }
Procedure Rectangle(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a + dx/2;
f := fx(c);
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления интеграла ');
readln(eps);
Number(a, b, eps, n);
Rectangle(a, b, n, j);
writeln('Значение интеграла равно ', j:6:t(eps));
writeln('С точностью до ', eps:1:t(eps))
end.
3.2. Формула трапеций
Рис. 49
Заменим данную кривую вписанной в нее ломаной, с вершинами в точках (xi, yi), где yi = f(xi) (i = 0, 1, 2, ..., n-1). Тогда криволинейная трапеция заменится фигурой, состоящей из трапеций (см. рис. 49). Будем по-прежнему считать, что промежуток [a, b] разбит на равные части, тогда площади этих трапеций будут равны:
Складывая полученные значения, приходим к приближенной формуле:
Эта приближенная формула называется формулой трапеций.
Оценка погрешности формулы трапеций определяется по следующему дополнительному члену:
Задание 4
Составьте процедуры вычисления интеграла по формуле трапеций и процедуру определения числа точек деления.
Составьте полностью программу и вычислите интегралы следующих функций на указанных промежутках с заданной точностью:
1) f(x) = 1/x на [1; 2] с точностью до 0.001;
2) f(x) = x/(x2
+ 1)2 на [-1; 2] с точностью до 0.01;
3) f(x) = 1/cos2 x на [0; Pi/3] с точностью до 0.0000001.
3.3. Параболическое интерполирование. Дробление промежутка интегрирования. Формула Симпсона
Параболическое интерполирование
Для приближенного вычисления интеграла функции
и тогда будет выполняться приближенное равенство
Такая замена тем более важна, когда подынтегральная функция представляет собой выражение, интеграл которой точно вычислить нельзя.
Геометрически это можно представить так, что криволинейная трапеция, находящаяся под кривой y = f(x) заменяется "параболой k-го порядка", поэтому такой процесс получил название параболическое интерполирование.
Существует много интерполяционных многочленов, о некоторых из них мы будем говорить на других занятиях, а на этом познакомимся с интерполяционным многочленом Лагранжа, который строится по следующей формуле:
При интегрировании получается линейное относительно значений
При k = 0, функция f(x) просто заменяется постоянной f(c0), где c0 - любая точка из промежутка [a, b], например средняя: c0 = (a + b)/2. Тогда приближенно
Геометрически это означает, что криволинейная фигура под кривой y = f(x) заменяется прямоугольником с высотой, равной средней ординате.
При k = 1 функция f(x) заменяется линейной функцией P1(x), которая имеет одинаковые с ней значения при x = c0 и x = c1. Если взять c0
= a, c1 = b, то
после преобразования, получим
Таким образом, здесь мы приближенно допускаем
Геометрически это будет представляться, как замена криволинейной фигуры трапецией: вместо кривой берется хорда, соединяющая ее концы.
Более интересный результат получается, если взять k = 2. Можно положить c0= a, c1 = (a + b)/2, c2 = b, тогда интерполяционный многочлен P2(x), будет иметь вид
После преобразований и интегрирования многочлена P2(x), приходим к приближенной формуле
В этом случае, площадь фигуры под кривой заменяется площадью фигуры, ограниченной параболой (с вертикальной осью), проходящей через крайние и среднюю точки кривой.
Увеличивая степень k интерполяционного многочлена, т.е. проводя параболу через все большее число точек данной кривой, можно добиться большей точности. Но на практике часто используют другой способ, основанный на сочетании параболического интерполирования с идеей дробления промежутка.
3.4. Дробление промежутка интегрирования
Для вычисления интеграла
[x0, x1], [x1, x2], ..., [xn-1, xn] (x0 = a, xn = b),
искомый интеграл представится в виде суммы
После этого, к каждому из этих промежутков применим параболическое интерполирование, т.е. станем вычислять перечисленные интегралы по одной из приближенных формул: прямоугольников, трапеций или по параболической формуле (7).
В первых двух случаях мы получим уже известные формулы прямоугольников и трапеций.
Применим теперь формулу (7) к каждому из интегралов, при этом положим, что
Получим
................................................
Складывая почленно эти равенства, получим формулу:
Эта формула называется формулой Симпсона. Ею пользуются для приближенного вычисления интегралов чаще, чем формулами прямоугольников и трапеций, так как она дает более точный результат.
Составим процедуру на языке Паскаль реализующую эту формулу.
{ Вычисление интеграла по формуле Симпсона }
Procedure Simpson(a, b : real; n : integer; var j : real);
var
dx, c, c1, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a;
c1 := a + dx/2;
f := fx(a) + fx(b) + 4*fx(c1);
for i := 1 to n - 1 do
begin
c := c + dx; c1 := c1 + dx;
f := f + 2*fx(c) + 4*fx(c1)
end;
j := (dx/6)* f
end;
Остается выяснить вопрос о числе точек деления в зависимости от заданной точности вычисления, другими словами, установить погрешность вычисления.
Если промежуток [a, b] разделен на n равных частей, то для формулы Симпсона дополнительный член имеет вид
Значит оценить точность вычисления можно по модулю этого дополнительного члена: |Rn|.
Однако, это связано с неприятностью находить 4-ю производную функции, что достаточно нелегкое дело. Есть и другой путь оценки точности вычисления интегралов и не только по формуле Симпсона, но и по формуле прямоугольника и трапеций.
Этот прием заключается в следующем. Искомый интеграл вычисляется дважды: при делении отрезка [a, b] на n частей и на 2n частей. Полученные значения интегралов In
и I2n сравниваются и первые совпадающие десятичные знаки считаются верными.
Покажем, как, используя такой метод оценить точность интеграла, вычисленного по формуле Симпсона.
Преобразуем дополнительный член Rn формулы Симпсона:
Пусть Rn
и R2n - погрешности интегрирования, тогда, учитывая предыдущую формулу можно составить пропорцию:
Понятно, что h2n
= hn/2. Тогда из пропорции получаем: Rn = 16 R2n. Если I - истинное значение интеграла, то I = In + Rn и I = I2n + R2n, откуда находим: In + 16 R2n
= I2n + R2n, т.е.
Надо помнить, что в процессе практических вычислений при последовательном удвоении числа отрезков разбиения начинает сильно прогрессировать удельный вес ошибки округления, значение которой с некоторого момента ставит предел достижимой точности.
Ниже, мы рассмотрим и другой способ оценки точности. А сейчас, основываясь на последней формуле, составим программу вычисления интеграла с применением формулы Симпсона.
{ Вычисление интеграла по формуле Симпсона }
Program Jntegral_Simpson;
uses WinCrt;
var
a, b, eps, j, j1 : real;
n : integer;
Function fx(x : real) : real;
begin
fx := exp(-x*x)
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Вычисление интеграла по формуле Симпсона }
Procedure Simpson(a, b : real; n : integer; var j : real);
var
dx, c, c1, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a;
c1 := a + dx/2;
f := fx(a) + fx(b) + 4*fx(c1);
for i := 1 to n - 1 do
begin
c := c + dx; c1 := c1 + dx;
f := f + 2*fx(c) + 4*fx(c1)
end;
j := (dx/6)* f
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления интеграла '); readln(eps);
n := 2;
Simpson(a, b, n, j); Simpson(a, b, 2*n, j1);
while abs(j - j1)/15 > eps do
begin
n := n + 2;
Simpson(a, b, n, j); Simpson(a, b, 2*n, j1)
end;
writeln('Значение интеграла равно ', j1:6:t(eps));
writeln('С точностью до ', eps:1:t(eps))
end.
3.5. Об оценке погрешности
При вычисление интегралов по формуле прямоугольников и по формуле трапеций оценка погрешности делалась по формуле с помощью второй производной, которая была ограничена на заданном промежутке своим наибольшим значением.
Так, в формуле трапеций, оценка погрешности использовалась для определения числа точек разбиения, которое определяло точность вычисления. Это делалось с помощью процедуры:
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps, max : real; var
n : integer);
var
d : real;
begin
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while (max*d)/(24*n*n) >= eps do n := n+1;
end;
Такого рода оценки погрешности называются гарантированными оценками погрешности.
На основании этой оценки можно гарантировать, что погрешность приближенного значения интеграла не превосходит определенной величины.
Существует второй способ оценки - отказ от получения строгой, гарантированной оценки погрешности и получение оценки погрешности лишь с определенной степенью достоверности. В частности при оценки погрешности вычисления интеграла по методу Симпсона погрешность оценивалась через разность результатов приближенного значения интеграла при различных значениях параметров, в частности, при n и 2n.
Необходимость применения метода Симпсона, как раз по причине метода оценки погрешности, примененной в этом методе, в следующих случаях.
1) Когда уже первая производная подынтегральной функции равна нулю.
2) Когда первая производная подынтегральной функции обращается в точках (пределах интегрирования или других) в бесконечность.
Примером может быть интеграл
Подынтегральная функция при
Эти два обстоятельства вызывают необходимость выбирать метод интегрирования в зависимости от заданной функции и поведения ее производных.
Задание 5
1. Используя оператор выбора Case ... of ..., составьте программу, с помощью которой пользователю можно было бы выбирать метод интегрирования, подобно программе выбора метода решения уравнений.
2. Вычислите по формуле Симпсона следующие интегралы:
1) полный эллиптический интеграл 2-го рода
2)
3.6. Вычисление интегралов методом Монте-Карло
Вычисление интегралов методом Монте-Карло часто применяется для двойных, тройных, вообще, кратных интегралов. Идея метода состоит в следующем.
Пусть задана некоторая подынтегральная функция F - непрерывная в области интегрирования Q. Выберем в этой области n случайных точек M, найдем значение заданной функции в некоторой "средней" точке области интегрирования. При достаточно большом n можно считать, что
Тогда, значение интеграла приблизительно равно
Применим этот метод для простейшего интеграла на промежутке [a, b], т.е. необходимо вычислить интеграл:
В этом случае в качестве объема области интегрирования D выступает длина отрезка [a, b], которая равна: D = b - a.
Пусть xi (i = 1, 2, ..., n) - случайные точки из промежутка [a, b], тогда значение функции f(x) в некоторой "средней" точке будет:
а значение интеграла станет равно
Для получение точек xi можно использовать уже известный способ нахождения случайных точек с помощью функции random,
x := random*(b - a) + a
Функция для вычисления интеграла получится такой:
{ Функция вычисления интеграла методом Монте-Карло }
Function I(n : longint; a, b : real) : real;
var
x, f : real;
k : longint;
begin
randomize;
f := 0;
for k := 1 to n do
begin
x := random*(b - a) + a;
f := f + fx(x)
end;
I := (b - a)*f/n
end;
Итак, с помощью генератора случайных чисел (random) вырабатывается случайное число из промежутка [a, b], находится значение функции в этой точке и суммируются ее абсолютные величины (учитывая, что функция может принимать и отрицательные значения) в переменную f. Приближенное значение интеграла вычисляется с помощью оператора I := (b - a)*f/n.
Теперь стоит очень болезненный вопрос оценки точности значения интеграла. Не вдаваясь в долгие математические рассуждения (о чем смотрите [19]), покажем неравенства, с помощью которых можно оценить точность вычисления интеграла.
Пусть I(f) - точное значение интеграла, а Sn(f) - его приближенное значение, тогда с вероятностью 0,997 и 0,99999 выполняются следующие соотношения
Здесь D(f) - дисперсия непрерывной случайной величины, которая вычисляется по формуле:
где сумма есть математическое ожидание случайной величины.
Приведем схему последовательного вычисления интеграла с заданной точностью eps. Последовательно, при n = 1, ... получаются случайные точки x и вычисляющая величины tn
, Sn , dn , пользуясь рекуррентными соотношениями
и величину
Начальные условия рекурсии n = 1,
Если оказалось, что
Изложенные математические соображения можно реализовать следующей процедурой:
{ Процедура оценки точности и вычисления интеграла }
Procedure Monte_Karlo(a, b, eps : real; var
s : real);
var
t, d, f, dd : real;
n : longint;
begin
n := 1;
t := I(n, a, b);
s := t;
d :=0;
repeat
n := n + 1;
t := t + I(n, a, b);
s := t/n;
d := d + (n/(n - 1))*sqr(I(n, a, b) - s);
dd := d/(n - 1)
until 5*sqrt(dd/n) < eps
end;
Полностью программа приводится ниже. Надо заметить, что для интегралов на промежутке программа работает медленно и лучше применять для их вычисления другие методы, а вот для кратных интегралов программа может быть очень полезной.
{Вычисление интеграла и оценка его точности методом Монте-Карло}
Program Integral_Monte_Karlo;
uses WinCrt;
var
a, b, s, eps : real;
{-----------------------------------------------------------------------------------------}
Function fx(x : real) : real; { Промежуток интегрир. [3, 4] }
begin { Интеграл хорошо вычисл. с точн. до 0.00001 }
fx := sin(0.2*x - 3)/(x*x + 1)
end;
{-----------------------------------------------------------------------------------------}
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления интеграла методом Монте-Карло }
Function I(n : longint; a, b : real) : real;
var
x, f : real;
k : longint;
begin
randomize;
f := 0;
for k := 1 to n do
begin
x := random*(b - a) + a;
f := f + fx(x)
end;
I := (b - a)*f/n
end;
{----------------------------------------------------------------------------------------}
{ Процедура вычисления точности и интеграла }
Procedure Monte_Karlo(a, b, eps : real; var s : real);
var
t, d, f, dd : real;
n : longint;
begin
n := 1;
t := I(n, a, b);
s := t;
d :=0;
repeat
n := n + 1;
t := t + I(n, a, b);
s := t/n;
d := d + (n/(n - 1))*sqr(I(n, a, b) - s);
dd := d/(n - 1)
until 5*sqrt(dd/n) < eps
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
randomize;
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления '); readln(eps);
Monte_Karlo(a, b, eps, s);
writeln('Значение интеграла равно ', s:6:t(eps));
write('С точностью до ', eps:1:t(eps));
writeln(' и вероятностью 0.99999')
end.
Задание 6
Используя второй способ оценки погрешности через разность результатов приближенного значения интеграла при различных значениях аргументов, изменить программу вычисления интеграла по методу Монте-Карло. И вычислить с ее помощью интегралы из задания 5.
3.7. Вычисление двойных интегралов методом Монте-Карло
Вначале сделаем отступление в курс математического анализа и познакомимся с понятием двойного интеграла. Как и прежде, в определении понятий математического анализа, мы будем следовать Г.М. Фихтенгольцу.
Возникновение двойного (определенного) интеграла связывают с задачей определения объема цилиндрического бруса, подобно тому, как определенный интеграл возник из задачи определения площади криволинейной трапеции.
Рассмотрим тело (V), которое сверху ограничено поверхностью z = f(x, y), с боков - цилиндрической поверхностью с образующими, параллельными оси z, наконец, снизу - плоской фигурой (P) на плоскости xy. Требуется найти объем V тела.
Для решения этой задачи мы прибегаем к обычному в интегральном исчислении приему, состоящему в разложении искомой величины на элементарные части, приближенному подсчету каждой части, суммированию и последующему предельному переходу. С этой целью разложим область (P) сетью кривых на части (P1), (P2), ..., (Pn) и рассмотрим ряд цилиндрических столбиков, которые имеют своими основаниями эти частичные области и в совокупности составляют данное тело.
Для подсчета объема отдельных столбиков возьмем произвольно в каждой фигуре (Pi) по точке (xi, yi). Если приближенно принять каждый столбик, за настоящий цилиндр с высотой, равной аппликате f(xi, yi), то объем отдельного столбика оказывается приближенно равным f(xi, yi)*Pi, где Pi означает площадь фигуры (Pi). В таком случае приближенное выражение объема всего тела будет
Для повышения точности этого равенства будем уменьшать размеры площадок (Pi), увеличивая их число. В пределе, при стремлении к нулю наибольшего из диаметров всех областей (Pi), это равенство делается точным, так что
и поставленная задача решена.
Предел этого вида и есть двойной интеграл от функции f(x, y) по области (P); он обозначается символом
Таким образом, двойной интеграл является прямым обобщением понятия простого определенного интеграла на случай функции двух переменных.
Пример 1. Вычислить интеграл, распространенный на прямоугольник
(P) = [3, 4; 1, 2]:
Решить этот интеграл методом Монте-Карло очень просто. Для этого достаточно изменить функцию:
Function fx(x, y : real) : real;
begin
fx := 5*x*x*y - 2*y*y*y
end;
Изменить функцию вычисления интеграла, куда добавить в качестве входных параметров две другие координаты прямоугольной области по оси OY и в результате вычисления умножать среднее значение функции, не на длину отрезка, а на площадь области, в данном случае, на площадь прямоугольника. Разумеется, задавать случайные значения для y из соответствующего промежутка [a1, b1].
Функция станет такой:
{ Функция вычисления интеграла методом Монте-Карло }
Function I(n : longint; a, b, a1, b1 : real) : real;
var
x, y, f : real;
k : longint;
begin
randomize;
f := 0;
for k := 1 to n do
begin
y := random*(b1 - a1) + a1;
x := random*(b - a) + a;
f := f + fx(x, y)
end;
I := (b - a)*(b1 - a1)*f/n
end;
Необходимо внести незначительные изменения в процедуру Monte_Karlo, а в основной программа, не забыть описать новые переменный и ввести их в программу.
Составьте программу самостоятельно и выполните ее.
Задание 7
Вычислить двойной интеграл
где (P) есть круг радиуса R с центром в начале координат.
Библиотека часто встречающихся процедур и функций
43. Процедура уточнения корня методом хорд.
{ Процедура уточнения корня методом хорд }
Procedure chord(a, b, eps, min : real; var
x : real);
var
x1 : real;
begin
x1 := a;
repeat
x := x1 - ((b - x1)*fx(x1))/(fx(b) - fx(x1));
x1 := x
until abs(fx(x))/min < eps
end;
44. Функция вычисления первой производной.
{ Вычисление 1-й производной и опред. точности ее вычислен.}
{ derivative - производная }
Function derivat1(x0, eps : real) : real;
var
dx, dy, dy2 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx/2) - fx(x0 - dx/2);
dy2 := fx(5*x0/4 + dx) - 2*fx(5*x0/4);
dy2 := dy2 + fx(5*x0/4 - dx)
until abs(dy2/(2*dx)) < eps;
derivat1 := dy/dx
end;
45. Процедура определения наименьшего значения первой производной.
{ Процедура определения наименьшего значения производной }
{ на заданном промежутке }
Procedure minimum(a, b, eps : real; var min : real);
var
d : real;
begin
a := a - eps;
b := b + eps;
repeat
a := a + eps;
b := b - eps;
min := abs(derivat1(a, eps));
d := abs(derivat1(b, eps));
if min > d then min := d
until min <> 0
end;
46. Функция вычисления порядка, в зависимости от задаваемой точности.
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
47. Процедура уточнения корня методом касательных.
{ Процедура уточнения корня методом касательных }
Procedure tangent(a, b, eps, min, dy : real; var x : real);
var
x1 : real;
begin
x1 := a;
repeat
x := x1 - fx(x1)/derivat1(x1);
x1 := x
until abs(fx(x))/min < eps
end;
48. Функция вычисления второй производной.
{ Функция вычисления второй производной }
Function
derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
49. Процедура вычисления корня уравнения (комбинированный метод).
{ Комбинированный метод }
Procedure Combination(a, b, eps : real; var x : real);
var
z : real;
begin
repeat
if fx(a)*derivat2(a, eps) > 0
then
begin
tangent(a, b, eps, z);
chord(b, a, x);
b := z; a := x
end
else
begin
tangent(b, a, eps, z);
chord(a, b, x);
b := x; a := z
end
until abs(z - x) < eps
end;
50. Процедура определения числа точек деления промежутка в зависимости от вычисляемой точности.
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps : real; var
n : integer);
var
dy2, d, c, dx : real;
begin
c := (a + b)/2;
dy2 := derivat2(c, eps);
if dy2 = 0
then
begin
c := a; dx := (b - a)/10;
while derivat2(c, eps) = 0 do c := c + dx;
dy2 := derivat2(c, eps)
end;
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while abs(dy2*d)/(24*n*n) >= eps do n := n+1;
end;
{--------------------------------------------------------------}
{ Вычисление интеграла методом прямоугольников }
Procedure Rectangle(a, b : real; n : integer; var j : real);
var
dx, c, f : real; i : integer;
begin
dx := (b - a)/n; c := a + dx/2;
f := fx(c);
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
51. Процедура вычисления интеграла по формуле прямоугольников.
{ Вычисление интеграла методом прямоугольников }
{ Rectangle - прямоугольник }
Procedure Rectangle(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a + dx/2;
f := fx(c);
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
52. Процедура вычисления интеграла по формуле трапеций.
{ Вычисл. интеграла по формуле трапеций. Trapezoid - трапеция }
Procedure Trapezoid(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n; c := a;
f := (fx(a)) + fx(b))/2;
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
53. Процедура вычисления интеграла по формуле Симпсона.
Procedure Simpson(a, b : real; n : integer; var j : real);
var
dx, c, c1, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a;
c1 := a + dx/2;
f := fx(a) + fx(b) + 4*fx(c1);
for i := 1 to n - 1 do
begin
c := c + dx; c1 := c1 + dx;
f := f + 2*fx(c) + 4*fx(c1)
end;
j := (dx/6)* f
end;
54. Процедуры вычисления интеграла по методу Монте-Карло.
Function I(n : longint; a, b : real) : real;
var
x, f : real; k : longint;
begin
randomize;
f := 0;
for k := 1 to n do
begin
x := random*(b - a) + a;
f := f + fx(x)
end;
I := (b - a)*f/n
end;
{----------------------------------------------------------------------------------------}
Procedure Monte_Karlo(a, b, eps : real; var s : real);
var
t, d, f, dd : real; n : longint;
begin
n := 1;
t := I(n, a, b);
s := t;
d :=0;
repeat
n := n + 1;
t := t + I(n, a, b);
s := t/n;
d := d + (n/(n - 1))*sqr(I(n, a, b) - s);
dd := d/(n - 1)
until 5*sqrt(dd/n) < eps
end;
Упражнения
184. Вычислить значение дифференциала функции
185.
186. Доказать, что функция y = ex sinx удовлетворяет соотношению
187. Решите уравнения, выбирая подходящий метод решения:
x3 - 9x + 2 = 0, xeч = 2, x = 0.538sinx + 1, aч = ax при a > 1, x2 arctgx = a, где
188. Решить уравнение
189. Пользуясь правилами прямоугольника, трапеции и правилом Симпсона, вычислить приближенно число
190. Вычислить
Найти модуль перехода от натурального логарифма к десятичному.
191. Вычислить по формуле Симпсона или методом Монте-Карло.
192. В следующих задачах при нахождении пределов интегрирования необходимо воспользоваться методами приближенного решения уравнений.
а) Найти площадь фигуры, ограниченной дугами парабол y = x3- 7 и
y = - 2x2+ 3x и осью ординат.
б) Найти площадь фигуры, ограниченной параболой y = x3 и прямой
y = 7(x + 1).
в) Найти площадь фигуры, ограниченной параболой y=16 - x3 и полукубической параболой y =
193. Вычислить двойной интеграл по методу Монте-Карло:
если область (A) ограничена двумя параболами: y = x2 и y2 = x.
Ответы
К
заданию 1
Program derivative2;
uses WinCrt;
var
x0, eps, dx, dy, dy3 : real;
Function fx(x : real) : real;
begin
fx := x*x*x*x
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления второй производной }
Function derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите точку, в которой находится ');
write('вторая производная '); readln(x0);
write(' Введите точность вычисления второй производной ');
readln(eps);
dy := derivat2(x0, eps);
write('Вторая производная функции в точке ', x0:6:t(eps));
writeln(' равна ', dy:6:t(eps));
writeln('с точностью до ', eps:1:t(eps))
end.
К заданию 3
{ Выбор метода решения уравнений: метод хорд, метод касательн. }
{ комбинированный метод }
Program Case_method;
uses WinCrt;
var
a, b, x, eps : real;
k : integer;
{----------------------------------------------------------------------------------------}
{ Заданная функция }
Function fx(x : real) : real;
begin
fx := x*sin(x) - 0.5
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Вычисление 1-й производной }
Function derivat1(x0, eps : real) : real;
var
dx, dy, dy2 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx/2) - fx(x0 - dx/2);
dy2 := fx(5*x0/4 + dx) - 2*fx(5*x0/4);
dy2 := dy2 + fx(5*x0/4 - dx)
until abs((dy2*dy2*fx(x0))/(2*dx)) < eps;
derivat1 := dy/dx
end;
{----------------------------------------------------------------------------------------}
{ Вычисление 2-й производной }
Function derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
{----------------------------------------------------------------------------------------}
{ Процедура определения наименьшего значения производной }
{ на заданном промежутке }
Procedure minimum(a, b, eps : real; var min : real);
var
d : real;
begin
a := a - eps;
b := b + eps;
repeat
a := a + eps;
b := b - eps;
min := abs(derivat1(a, eps));
d := abs(derivat1(b, eps));
if min > d then min := d
until min <> 0
end;
{----------------------------------------------------------------------------------------}
{ Процедура уточнения корня методом хорд }
Procedure chord1(a, b, eps : real; var
x : real);
var
x1, min : real;
begin
minimum(a, b, eps, min);
x1 := a;
repeat
x := x1 - ((b - x1)*fx(x1))/(fx(b) - fx(x1));
x1 := x
until abs(fx(x))/min < eps
end;
{----------------------------------------------------------------------------------------}
{ Процедура уточнения корня методом касательных }
Procedure tangent1(a, b, eps, min : real; var x : real);
var
x1 : real;
begin
x1 := a;
repeat
x := x1 - fx(x1)/derivat1(x1, eps);
x1 := x
until abs(fx(x))/min < eps
end;
{----------------------------------------------------------------------------------------}
Procedure Tangent2(a, b, eps : real; var x : real);
var
min : real;
begin
minimum(a, b, eps, min);
if fx(a)*derivat2(a, eps) > 0
then tangent1(a, b, eps, min, x)
else tangent1(b, a, eps, min, x)
end;
{----------------------------------------------------------------------------------------}
{ Процедура уточнения корня методом хорд }
Procedure chord(a, b : real; var
x : real);
begin
x := a - ((b - a)*fx(a))/(fx(b) - fx(a))
end;
{----------------------------------------------------------------------------------------}
{ Процедура уточнения корня методом касательных }
Procedure tangent(a, b, eps : real; var
z : real);
begin
z := a - fx(a)/derivat1(a, eps)
end;
{----------------------------------------------------------------------------------------}
{ Комбинированный метод }
Procedure Combination(a, b, eps : real; var x : real);
var
z : real;
begin
repeat
if fx(a)*derivat2(a, eps) > 0
then
begin
tangent(a, b, eps, z);
chord(b, a, x);
b := z; a := x
end
else
begin
tangent(b, a, eps, z);
chord(a, b, x);
b := x; a := z
end
until abs(z - x) < eps
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите левый конец промежутка a = '); readln(a);
write('Введите правый конец промежутка b = '); readln(b);
write('Введите точность вычисления корня eps = ');
readln(eps);
writeln('Для выбора метода решения, введите его номер ');
writeln('1 - метод хорд, 2 - метод касательных');
writeln(' 3 - комбинированный метод'); readln(k);
Case k of
1 : chord1(a, b, eps, x);
2 : Tangent2(a, b, eps, x);
3 : Combination(a, b, eps, x)
end;
writeln('Корень уравнения равен x = ', x:6:t(eps));
writeln('С точностью до eps = ', eps:1:t(eps))
end.
К заданию 4
{ Вычисление интеграла по формуле трапеций }
Program integral_trapezoid;
uses WinCrt;
var
a, b, eps, j : real;
n : integer;
Function fx(x : real) : real;
begin
fx := 1/x
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления второй производной }
Function derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
{----------------------------------------------------------------------------------------}
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps : real; var
n : integer);
var
dy2, d, c, dx : real;
begin
c := (a + b)/2;
dy2 := derivat2(c, eps);
dx := (b - a)/10;
if dy2 = 0
then
begin
c := a;
while derivat2(c, eps) = 0 do c := c + dx;
dy2 := derivat2(c, eps)
end;
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while abs(dy2*d)/(12*n*n) >= eps do n := n+1;
end;
{----------------------------------------------------------------------------------------}
{ Вычисл. интеграла по формуле трапеций. Trapezoid - трапеция }
Procedure Trapezoid(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a;
f := (fx(a) + fx(b))/2;
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
{----------------------------------------------------------------------------------------}
begin
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления интеграла ');
readln(eps);
Number(a, b, eps, n);
Trapezoid(a, b, n, j);
writeln('Значение интеграла равно ', j:6:t(eps));
writeln('С точностью до ', eps:1:t(eps))
end.
2-й способ
{ Вычисление интеграла по формуле трапеций }
Program integral_trapezoid;
uses WinCrt;
var
a, b, eps, max, j : real;
n : integer;
Function fx(x : real) : real;
begin
fx := 1/x
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол-во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления второй производной }
Function derivat2(x0, eps : real) : real;
var
dx, dy, dy3 : real;
begin
dx := 1;
repeat
dx := dx/2;
dy := fx(x0 + dx) - 2*fx(x0) + fx(x0 - dx);
dy3 := fx(5*x0/4 + 2*dx) - 2*fx(5*x0/4 + dx);
dy3 := dy3 - fx(5*x0/4 - 2*dx) + 2*fx(5*x0/4 - dx)
until abs(dy3/(6*dx)) < eps;
derivat2 := dy/(dx*dx)
end;
{----------------------------------------------------------------------------------------}
{ Определение наибольшего значения второй производной }
Procedure Maximum(a, b, eps : real; var
max : real);
var
dx, x : real;
begin
dx := 0.1; x := a;
max := abs(derivat2(x, eps));
while x<= b do
begin
x := x + dx;
if max < abs(derivat2(x, eps))
then max := abs(derivat2(x, eps))
end
end;
{----------------------------------------------------------------------------------------}
{ Процедура определения числа точек деления промежутка интегр. }
Procedure Number(a, b, eps, max : real; var
n : integer);
var
d : real;
begin
n := 1;
d := abs((b - a)*(b - a)*(b - a));
while (max*d)/(24*n*n) >= eps do n := n+1;
end;
{----------------------------------------------------------------------------------------}
{ Вычисл. интеграла по формуле трапеций. Trapezoid - трапеция }
Procedure Trapezoid(a, b : real; n : integer; var j : real);
var
dx, c, f : real;
i : integer;
begin
dx := (b - a)/n;
c := a;
f := (fx(a) + fx(b))/2;
for i := 1 to n - 1 do
begin
c := c + dx;
f := f + fx(c)
end;
j := dx * f
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления интеграла '); readln(eps);
Maximum(a, b, eps, max); Number(a, b, eps, max, n);
Trapezoid(a, b, n, j);
writeln('Значение интеграла равно ', j:6:t(eps));
writeln('С точностью до ', eps:1:t(eps))
end.
К заданию 6
{ Вычисление интеграла методом Монте-Карло }
Program integral_Monte_Karlo;
uses WinCrt;
var
a, b, eps : real;
n : longint;
{----------------------------------------------------------------------------------------}
Function fx(x : real) : real; { Промежуток интегрир. [3, 4] }
begin { Интеграл хорошо вычисл. с точн. до 0.00001 }
fx := sin(0.2*x - 3)/(x*x + 1)
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления порядка - кол- во знаков после запятой }
Function t(eps : real) : integer;
var
k : integer;
begin
k := -1;
repeat
eps := eps*10;
k := k + 1
until eps > 1;
t := k
end;
{----------------------------------------------------------------------------------------}
{ Функция вычисления интеграла методом Монте-Карло }
Function I(n : longint; a, b : real) : real;
var
x, f : real;
k : longint;
begin
randomize;
f := 0;
for k := 1 to n do
begin
x := random*(b - a) + a;
f := f + fx(x)
end;
I := (b - a)*f/n
end;
{----------------------------------------------------------------------------------------}
{ Основная программа }
begin
write('Введите нижний предел интегрирования '); readln(a);
write('Введите верхний предел интегрирования '); readln(b);
write('Введите точность вычисления '); readln(eps);
n := 10;
repeat
n := 2*n
until abs(I(n, a, b) - I(2*n, a, b)) <= eps;
writeln('Значение интеграла равно ', I(2*n, a, b):6:t(eps));
writeln('С точностью до ', eps:1:t(eps));
end.
