6.3 Метод Ромберга

Можно вычислить приближение Симпсона по особым линейным комбинациям формул трапеций.

а b a b

T — площадь под T — площадь под

2⁰ = 1 трапецией 2¹ = 2 трапециями

T =

a b a b

T — площадь под T — площадь под

2² = 4 трапециями 2³ = 8 трапециями

Пусть промежуток интегрирования разбит на 2^к равных частей и для этого разбиение по формуле трапеций получено значение Т₀^(к). Значение Т₀^(к) совпадают со значениями вычисляемого интеграла, если интегрируемая функция линейна, то есть f (x) является многочленом 1-й степени. По формуле:

T_m^(k) = , называемой формулой Ромберга, построим следующую Т - схему:

T₀⁽⁰⁾

T₁⁽⁰⁾

T₀⁽¹⁾ T₂⁽⁰⁾

T₁⁽¹⁾ T₃⁽⁰⁾

T₀⁽²⁾ T₂⁽¹⁾ T₄⁽⁰⁾

T₁⁽²⁾ T₃⁽¹⁾

T₀⁽³⁾ T₂⁽²⁾

T₁⁽³⁾

T₀⁽⁴⁾

Для интегрируемых по Риману функций все столбцы и все наклонные строки Т — схемы сходятся к искомому значению интеграла. Выпишем явные формулы для следующего фрагмента Т — схемы:

T₀⁽⁰⁾

T₁⁽⁰⁾

T₀⁽¹⁾ T₃⁽⁰⁾

T₂⁽¹⁾

T₀⁽²⁾

1. k = 0; 2^k = 2⁰ = 1 часть, то есть (b - a). T₀⁽⁰⁾ = [f (x₀) + f (x₄)], так как при k = 2; 2^k = 2² = 4 части x_i = a + i * , i = 0, 1, 2, 3, 4.

2. T₀⁽¹⁾ = * [f (x₀) + f (x₂)] + * [f (x₂) + f (x₄)] = [f (x₀) + 2f (x₂) + f (x₄)]k = 1 → 2^k = 2¹ = 2 части

3. T₁⁽⁰⁾ = = = =

=

= .

4. T₀⁽²⁾ =

5. T₁⁽¹⁾ =

= = =

формула Симпсона (парабол)

6. T₂⁽⁰⁾ =

=

= -

формула Ньютона — Котеса

Задача: Известно, что

π

Определить приближённо число π, проведя для указанного интеграла 4 итерации по Ромбергу.

Дополнение 1:

3[ ] = 2

6

=

.

Так как и симметричен, то:

; ,

.

Дополнение 2. Изложим последовательность действий при решении задач методом Ромберга:

1. В качестве первого приближения вычисляется значение площади трапеции, основания которой проведены через границы промежутка интегрирования. Полученное значение заносится в матрицу приближений как:

.

2. Затем запускается цикл. На каждом его круге шаг уменьшается вдвое. На первом обороте цикла вычисляется приближение по методу трапеций из того условия, что интервалов интегрирования уже два. Полученное более точное, приближение заносится в первую строку нулевого столбца. На втором круге шаг уменьшается в четыре раза относительно первоначального, и соответственно трапеций будет просуммировано четыре. Новое приближение занесено во вторую строку. И так далее, до тех пор, пока не будет выполнено условие остановки цикла. В качестве условия остановки цикла в методе Ромберга используют критерий разности двух последних приближений. В этом случае, если разность по модулю окажется меньше заданной точности, то цикл остановится и в качестве ответа будет выдано последнее приближение.

Sub Rombeg() 'начало программы

Dim Result As Double 'задание типа параметра результата

Dim T() As Double 'задание типа параметра

Dim D As Double 'задание типа параметра длины отрезка

Dim S As Double 'задание типа параметра суммы

Dim H As Double 'задание типа параметра расстояния

Dim X As Double 'задание типа параметра х

Dim M As Double 'задание типа параметра

Dim N As Integer 'задание типа параметра

Dim I As Integer 'задание типа параметра переменной в цикле

Dim J As Integer 'задание типа параметра переменной в цикле

A = InputBox("введите нижний предел") 'диалог запроса нижнего предела

B = InputBox("введите верхний предел") 'диалог запроса верхнего предела

K = InputBox(" введите число разбиений") 'диалог запроса числа разбиений

ReDim T(1 To K + 1) 'объявление массива

D = B - A 'вычисление длины отрезка

T(1) = f(A) + f(B) / 2 'вычисление Т(1)

N = 1 'присвоение N значение 1

I = 1 'присвоение I значение 1

Do 'начало цикла

S = 0 'присвоение сумме значение 0

N = 2 * N 'увеличение N в двое

H = D / N 'вычисление расстояния

J = 1 'присвоение J значение 1

Do 'начало цикла

X = A + J * H 'вычисление х координаты

S = S + f(X) 'вычисление суммы

J = J + 2 'увеличение переменной на 2

Loop Until Not J <= N 'конец цикла с условием

T(I + 1) = (2 * S / N + T(I)) / 2 'вычисление последующего Т

M = 1 'присвоение M значение 1

J = I 'присвоение J значение I

Do 'начало цикла

M = 4 * M 'увеличение M в 4 раза

T(J) = T(J + 1) + (T(J + 1) - T(J)) / (M - 1) 'вычисление Т(J)

J = J - 1 'уменьшение J на 1

Loop Until Not J >= 1 'конец цикла с условием

I = I + 1 'увеличение I на 1

Loop Until Not I <= K 'конец цикла с условием

Result = T(1) * D 'вычисление результата

c = MsgBox(Result, 0, "Результат") 'выведение сообщения на экран

End Sub 'конец программы

Function f(X) 'объявление ф-ии

f = 4 / (1 + X * X) 'ф-ия задана уравнением

End Function 'конец объявления