Правило Рунге оценки погрешности

Приведем сводку квадратурных формул с остаточными членами [7].

Формула трапеций:

Формула Симпсона:

Рассмотрим на примере общей формулы Симпсона (6.16) правило Рунге для оценки погрешности. Пусть подынтегральная функция четырежды непрерывно дифференцируема. Тогда формула остаточного члена имеет вид:

(6.19)

где ξ — некоторое число из отрезка [ a, b ]. Предположим, что производная y ^IV(x) изменяется на этом отрезке медленно, и приближенно можно записать остаточный член в виде где M — постоянная. Пусть S_h и S _{2 h} соответственно значения интеграла, полученные по общей формуле Симпсона с шагом h и 2 h. Тогда справедливы соотношения

Отсюда получим формулу для оценки погрешности

. (6.20)

Уточненное значение интеграла по формуле Симпсона вычисляется с учетом поправки

(6.21)

Для формулы трапеций также можно применить правило Рунге. Так как формула остаточного члена общей формулы трапеций может быть представлена в виде , то справедливы соотношения

Отсюда получим формулу для оценки погрешности

. (6.22)

Теперь для интеграла можно записать по правилу Рунге формулу трапеций с поправкой

(6.23)

Замечание. Здесь необходимо подчеркнуть, что описанное правило Рунге применимо только тогда, когда выполняются указанные выше условия для функции y (x) — существование производной соответствующего порядка и её ограниченность, точнее говоря, возможность приближенного представления погрешности в виде для формулы Симпсона ( для формулы трапеций), где M — постоянная. Погрешность представления остаточного члена в указанном виде здесь считается достаточно малой. Эти условия для конкретной функции могут не выполняться, тогда правило Рунге не будет гарантировать приемлемый результат.

Пример 6.4. Вычислить интеграл по правилу Рунге, используя разбиение отрезка на n = 5 и n = 10 частей, используя формулу трапеций (6.14).

Решение в Excel. Используем результаты вычислений в примере 6.1. Для вычисления интеграла по формуле трапеций с разбиением отрезка [1, 2] на n = 5 частей получим таблицу 6.5.

Таблица 6.5

По формуле (6.23) получим

Найдем относительную погрешность этого значения

Применение правила Рунге позволило уменьшить относительную погрешность со значения 0, 001 до значения 0, 0000036.

Составим функции для вычисления интегралов с автоматическим выбором шага. С помощью меню «Сервис — Макроc — Редактор Visual Basic» откроем окно редактора, выполним команду «Insert — Module» и введем программы

Function f(x): f = 1 / x: End Function

Function Int_tr(a, b, n)

s = 0: h = (b - a) / n: x = a

For i = 1 To n - 1: x = x + h: s = s + f(x): Next i

Int_tr = h * (s + (f(a) + f(b)) / 2)

End Function

Function Int_tr_eps(a, b, eps)

n = 2: sh = Int_tr(a, b, n)

1 s2h = sh: n = 2 * n: sh = Int_tr(a, b, n): serr = (sh - s2h) / 3

If Abs(serr) > eps Then GoTo 1

Int_tr_eps = sh + serr

End Function

Введем в ячейку D 15 формулу =Int_tr_eps(1; 2; 0, 00001), получим значение 0, 693147181322587, которое имеет 8 верных знаков (точное значение интеграла 0, 693147180559945).

Приведем аналогичные функции для вычисления интеграла методом Симпсона с автоматическим выбором шага:

Option Explicit

Function f(ByVal x) As Double

f = 1 / x

End Function

Function Int_Simpson(ByVal a, ByVal b, ByVal n) As Double

Dim s, x, h, h2 As Double

Dim i As Integer

s = 0: h = (b - a) / n: h2 = h * 2: x = a

For i = 0 To n - 2 Step 2

s = s + f(x) + 4 * f(x + h) + f(x + h2): x = x + h2: Next i

Int_Simpson = h * s / 3: End Function

Function Int_Simpson_eps(ByVal a, ByVal b, ByVal eps) As Double

Dim sh, s2h, serr As Double

Dim n As Integer

n = 4: sh = Int_Simpson(a, b, n)

1 s2h = sh: n = 2 * n: sh = Int_Simpson(a, b, n): serr = (sh - s2h) / 15

If Abs(serr) > eps Then GoTo 1

Int_Simpson_eps = sh + serr: End Function

Введем в ячейку D 16 формулу =Int_Simpson_eps(1; 2; 0, 000000000001), получим результат, который совпадает с точным.

Замечание. Чтобы добиться максимально возможной точности, надо в функциях использовать переменные двойной точности.

Приведем программы на языке C ++ для вычисления интеграла методом Симпсона:

#include < iostream.h>

#include < math.h>

double f(double x);

typedef double (*PF)(double);

double Simpson(PF f, double a, double b, const int n);

double Simpson_eps(PF f, double a, double b, double eps);

int main(){

double a, b, y, eps; PF pf; int n;

cout < < " \n a = "; cin > > a; cout < < " \n b = "; cin > > b;

cout < < " \n n = "; cin > > n; cout < < " \n eps = "; cin > > eps;

pf = f;

y = Simpson(pf, a, b, n);

cout < < " \n n = " < < n< < " integral = " < < y;

y = Simpson_eps(pf, a, b, eps);

cout < < " \n eps = " < < eps< < " integral = " < < y;

cout < < " \n Press any key & Enter "; cin > > a;

return 0;

}

double f(double x){

double r;

r = sin(x);

return r;

}

double Simpson(PF f, double a, double b, const int n){

double s, h, h2, x; int i, n2;

s = 0; h = (b-a)/n; h2 = h*2; x = a; n2=n/2;

for (i = 0; i< n2; i++){

s = s+f(x)+4*f(x+h)+f(x+h2); x = x+h2; }

s = h*s/3;

return s;

}

double Simpson_eps(PF f, double a, double b, double eps){

typedef double (*PY)(double); PY py;

double s, sh, s2h, serr; int n;

n = 4; py = f; sh = Simpson(py, a, b, n);

do{ s2h = sh; n = 2*n; sh = Simpson(py, a, b, n);

serr = (sh-s2h)/15; }while(fabs(serr) > eps);

s = sh + serr;

return s;

}

Программа Simpson(pf, a, b, n) вычисляет интеграл по формуле Симпсона применяя деление отрезка интегрирования на n частей, а программа
Simpson_eps(pf, a, b, eps) применяет правило Рунге для метода Симпсона и вычисляет интеграл с точностью eps. В данном случае подынтегральная функция равна sin x.

Приведем результат счета: