Методы парабол (Симпсона) и более высоких степеней (Ньютона - Котеса)

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
Строка 103: Строка 103:
::<tex>\int_a^b{f(x)dx} = \sum_{i=0}^K{\int_{a_i}^{b_i}{f(x)dx}},</tex>
::<tex>\int_a^b{f(x)dx} = \sum_{i=0}^K{\int_{a_i}^{b_i}{f(x)dx}},</tex>
-
а на каждом из отрезков приближенное значение интеграла вычисляется при помощи квадратурных формул.
+
а на каждом из ''частичных отрезков'' <tex>[a_i,b_i]</tex> приближенное значение интеграла вычисляется при помощи квадратурных формул.
=== Примеры квадратурных формул ===
=== Примеры квадратурных формул ===
Строка 148: Строка 148:
== Рекомендации программисту ==
== Рекомендации программисту ==
 +
 +
=== Автоматический выбор шага интегрирования с помощью апостериорной оценки погрешности методом Рунге ===
 +
 +
Величина погрешности численного интегрирования зависит как от шага сетки <tex>h</tex>, так и от гладкости подынтегральной функции <tex>f(x)</tex>. В величину погрешности, кроме <tex>h</tex> входит также величина <tex>f^{(p)}(\xi)</tex>, которая может сильно меняться на отрезке <tex>[a,b]</tex> и заранее неизвестна. Для уменьшения величины погрешности, можно измельчить сетку на заданном отрезке. Но при этом необходимо апостериорно оценивать погрешность. Такую оценку погрешности можно осуществить ''методом Рунге''. Рассмотрим применение какой-то квадратурной формулы на частичном отрезке <tex>[x_{i-1},x_i]</tex>. Обозначим за <tex>I</tex> значение интеграла на всем отреpке, за <tex>I_i</tex> значение интеграла на <tex>i</tex>-ом частичном отрезке, за <tex>I_h</tex> приближенное значение интеграла на всем отрезке, полученное при помощи заданной квадратурной формулы и равномерном сетки с шагом <tex>h</tex>, а за <tex>I_{h,i}</tex> приближенное значение интеграла на <tex>i</tex>-ом частичном отрезке. Пусть данная квадратурная формула на данном частичном отрезке имеет порядок точности <tex>m</tex>, то есть
 +
 +
::<tex>I_i-I_{h,i} \approx ch^m,</tex>
 +
 +
где с <tex>c</tex> - некоторая константа. Тогда
 +
 +
::<tex>I_i-I_{ \frac{h}{2},i} \approx c( \frac{h}{2} )^m,</tex>
 +
 +
откуда получаем
 +
 +
::<tex>I_i-I_{h,i} \approx 2^m(I_i-I_{ \frac{h}{2},i})</tex>
 +
 +
{{ eqno | 8 }}
 +
::<tex>I_i-I_{ \frac{h}{2},i} \approx \frac {I_{ \frac{h}{2},i} - I_{h,i}} {2^m-1}.</tex>
 +
 +
Пусть используется составная квадратурная формула
 +
 +
::<tex>I \approx I_h=\sum_{i=1}^N{I_{h,i}},</tex>
 +
 +
причем на всех частичных отрезках используются квадратурные формулы с одним и тем же порядком точности (или же,в частности, одна и та же формула). Проведем на каждом частичном отрезке <tex>[x_{i-1},x_i]</tex> вычисления дважды - один раз с шагом <tex>h_i</tex>, второй раз с шагом <tex>\frac{h_i}{2}</tex> и апостериорно оценим погрешность по правилу Рунге {{ eqref | 8 }}. Если для заданного <tex>\varepsilon >0</tex> при <tex>i=1,2, \ldots N</tex> будут выполняться неравенства
 +
 +
{{eqno | 9 }}
 +
::<tex>|I_i-I_{ \frac{h}{2},i}| \approx \frac {|I_{ \frac{h}{2},i} - I_{h,i}|} {2^m-1} \le \frac {\varepsilon h_i} {b-a},</tex>
 +
 +
то получим
 +
 +
::<tex>|I-I_{\frac{h}{2}}| \le \frac {\varepsilon}{b-a} \sum_{i=1}^N{h_i} \le \varepsilon,</tex>
 +
 +
то есть будет достигнута заданная точность <tex>\varepsilon</tex>.
 +
 +
Если же на каком-то из частичных отрезков оценка {{ eqref | 9 }} не будет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида {{ eqref | 9 }}. Заметим, что для некоторой функции <tex>f(x)</tex> такое измельчение может продолжаться слишком долго. Поэтому в соответствующей программе следует предусмотреть ограничение сверху на число измельчений.
 +
 +
Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции <tex>f(x)</tex> и с мелким шагом - на участках быстрого изменения <tex>f(x)</tex>. Это позволяет при заданной точности <tex>\epsilon</tex> уменьшить количество вычислений значений <tex>f(x)</tex> по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем, что для нахождения сумм <tex>I_{ \frac{h}{2},i}</tex> не надо пересчитывать значения <tex>f(x)</tex> во всех узлах, достаточно вычислять <tex>f(x)</tex> только в новых узлах.
== Заключение ==
== Заключение ==

Версия 17:09, 18 октября 2008

Содержание

Введение

Постановка математической задачи

Задача численного интегрирования состоит в нахождении приближенного значения интеграла

( 1 )
J[f]=\int_a^b{f(x)dx},

где f(x) - заданная и интегрируемая на отрезке [a,b] функция. На отрезке вводится сетка \omega=\{x_i:x_0=a<x_1<\ldots<x_i<\ldots<x_N=b\} и в качестве приближенного значения интеграла рассматривается число

( 2 )
J_N[f]=\sum_{i=0}^N {c_i f(x_i)},

где f(x_i) - значения функции f(x) в узлах x=x_i , где c_i - весовые множители, зависящие только от узлов, но не зависящие от выбора f(x). Формула (2) называется квадратурной формулой. Задача численного интегрирования при помощи квадратур состоит в отыскании таких узлов \{x_i\} и таких весов \{c_i\}, чтобы погрешность квадратурной формулы

D[f]=\sum_{i=0}^N{c_i f(x_i)} - \int_a^b{f(x)dx} = J_N[f] - J[f]

была минимальной по модулю для функции из заданного класса (величина D[f] зависит от гладкости f(x)). Погрешность зависит как от расположения узлов, так и от выбора весовых коэффициентов. Введем на [a,b] равномерную сетку с шагом h, т.е. множество точек \omega_h=\{x_i=a+ih, i=0,1,\ldots,N,hN=b-a}, и представим интеграл (1) в виде суммы интегралов по частичным отрезкам:

( 3 )
\int_a^b{f(x)dx}=\sum_{i=1}^N{\int_{x_{i-1}}^{x_i}{f(x)dx}}.

Для построения формулы численного интегрирования на всм отрезке [a,b] достаточно построить квадратурную формулу для интеграла

( 4 )
\int_{x_{i-1}}^{x_i}{f(x)dx}

на частичном отрезке [x_{i-1},x_i] и воспользоваться свойством (3).

Построение квадратурных формул

В силу вышеизложенного выше, вычисление приближенного значения интеграла производится при помощи квадратурной формулы

\int_a^b{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.

Данную формулу при помощи замены x=a+(b-a)s, f(x)=f(a+(b-a)s) можно привести к стандартному виду

( 5 )
\int_0^1{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.

В общем случае узлы и веса неизвестны и подлежат определению.

Рассмотрим случай, когда узлы заданы и требуется найти веса квадратурной формулы \{c_k\}. Будем пользоваться требованием: формула (5) должна быть точной для любого полинома P_r(s) степени r\le m. Для того чтобы полином степени r удовлетворял данному требованию, достаточно потребовать, чтобы квадратурная формула была точной для любого одночлена s^\sigma степени \sigma (\sigma=0,1,\ldots,r). Учитывая, что J[s^\sigma]=\frac{1}{\sigma+1}, получаем m+1 уравнение

c_0+c_1+\ldots+c_m=1,
c_0 s_0+c_1 s_1+\ldots+c_m s_m=\frac{1}{2},


.........................................


c_0 s_0^\sigma + c_1 s_1^\sigma + \ldots + c_m s_m^\sigma=\frac{1}{\sigma+1},


.........................................


c_0 s_0^m + c_1 s_1^m + \ldots + c_m s_m^m=\frac{1}{m+1}.

Эта система имеет единственное решение, так как ее определителем является определитель Вандермонда, отличный от нуля если нет совпадающих узлов, s_0<s_q<\ldots<s_m.

Так, полагая m=2,s_0=0,s_1=\frac{1}{2},s_2=1, имеем систему c_0+c_1+c_2=1, \frac{c_1}{2}+c_2=\frac{1}{2}, \frac{c_1}{4}+c_2=\frac{1}{3}, решением которой являются веса формулы Симпсона: c_0=c_2=\frac{1}{6},c_1=\frac{4}{6}. Таким образом, формула Симпсона является точной для полинома второй степени. Однако, в силу симметрии, она является точной и для всех полиномов третьей степени:

P_3(s)=1+\alpha_1 (s-\frac{1}{2})+\alpha_2 (s-\frac{1}{2})^2 +\alpha_3 (s-\frac{1}{2})^3,

так как она точна для f(s)=(s-\frac{1}{2})^3:

J_2[(s-\frac{1}{2})^3]=\frac{1}{6}((-\frac{1}{2})^3+4\cdot 0+(\frac{1}{2})^3)=0,
L[(s-\frac{1}{2})^3]=\int_0^1{(s-\frac{1}{2})^3ds}=0.

Формулы треугольника и трапеции точны для линейной функции,т.е. для полинома первой степени, в чем легко убедиться непосредственно. В общем случае в качестве P_m(s) можно выбрать интерполяционный полином Лагранжа

P_m(s)=\sum_{k=0}^m{l_k^{(m)}(s)f(s_k)},

где l_k^{(m)}(s) - интерполяционный коэффициент Лагранжа. Из равенства

L[P_m]=\int_0^1{P_m(s)ds}=\sum_{k=0}^m{f(s_k)\int_0^1{l_k^{(m)}(s)ds}}=\sum_{k=0}^m{c_k f(s_k)}

видно, что формула (5) верна для полинома степени m, если весовые коэффициенты c_k определяются по формуле

( 6 )
c_k=\int_0^1{l_k^{(m)}(s)ds}.

Формулы такого типа называют квадратурными формулами Котеса.

Изложение метода

Применение квадратурных формул

Вернемся к рассмотрению интеграла (1). Как было показано выше, данный интеграл заменой сводится к интегралу на единичном отрезке [0,1], а следовательно легко обобщить формулы для приближенного вычисления интеграла на единичном отрезке на произвольный. Применим аппарат квадратурных формул. Пусть задано равномерное разбиение отрезка с шагом h: \omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}, где обозначим f_i=f(a+ih). Пусть также выбрана некоторая квадратурная формула Ньютона-Котеса (то есть выбрана степень полинома n, а следовательно каждый полином строится по n+1 точке сетки). Мы также считаем, что данный набор из N+1 точки можно разбить на поднаборы по m точек с совпадающими крайним, то есть

( 7 )
N=n*k, где k \in \mathbb N.

Тогда, суммируя значения квадратур на каждом поднаборе, получим приближенное значение искомого интеграла. Если обозначить весовые множители c_i,i=0 \ldots n, то приближенное значение интеграла можно записать в виде двойной суммы

\int_a^b{f(x)dx \approx \sum_{j=0}^{k-1}{\sum_{i=0}^n{c_i f_{(i+j*n)}}}

Данный алгоритм естественным образом обобщается на случай, когда [a,b]=\bigcup_{i=0}^K{[a_i,b_i]}, где a_0=a, b_K=b, b_{i+1}=a_i;i=0\ldots K-1, а на каждом отрезке [a_i,b_i] задана равномерная сетка. Тогда искомый интеграл равен

\int_a^b{f(x)dx} = \sum_{i=0}^K{\int_{a_i}^{b_i}{f(x)dx}},

а на каждом из частичных отрезков [a_i,b_i] приближенное значение интеграла вычисляется при помощи квадратурных формул.

Примеры квадратурных формул

Приведем примеры квадратурных формул Котеса на равномерной сетке с шагом h: \omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}, где обозначим f_i=f(a+ih), \xi \in [a,b]:

  • для 3 точек(метод Симпсона), m=2:
J_2[f]=\frac{1}{3}h(f_0+4f_1+f_2),
D[f]=\frac{1}{90} h^5 f^{(4)}(\xi),
  • для 4 точек, m=3:
J_3[f]=\frac{3}{8}h(f_0+3f_1+3f_2+f_3),
D[f]=\frac{3}{80} h^5 f^{(4)}(\xi),
  • для 5 точек, m=4:
J_4[f]=\frac{2}{45}h(7f_0+32f_1+12f_2+32f_3+7f_4),
D[f]=\frac{8}{945} h^7 f^{(6)}(\xi),
  • для 6 точек, m=5:
J_5[f]=\frac{5}{288}h(19f_0+75f_1+50f_2+50f_3+75f_4+19f_5),
D[f]=\frac{275}{12096} h^7 f^{(6)}(\xi),
  • для 7 точек, m=6:
J_6[f]=\frac{1}{140}h(41f_0+216f_1+27f_2+272f_3+27f_4+216f_5+41f_6),
D[f]=\frac{9}{1400} h^9 f^{(8)}(\xi),
  • для 8 точек, m=7:
J_7[f]=\frac{7}{17280}h(751f_0+3577f_1+1323f_2+2989f_3+2989f_4+1323f_5+3577f_6+751f_7),
D[f]=\frac{8183}{518400} h^9 f^{(8)}(\xi),
  • для 9 точек, m=8:
J_8[f]=\frac{4}{14175}h(989f_0+5888f_1-928f_2+10496f_3-4540f_4+10496f_5-928f_6+5888f_7+989f_8),
D[f]=\frac{2368}{467775} h^{11} f^{(10)}(\xi),
  • для 10 точек, m=9:
J_9[f]=\frac{9}{89600}h(2857f_0+15741f_1+1080f_2+19344f_3+5778f_4+5778f_5+19344f_6+1080f_7+15741f_8+2857f_9).
D[f]=\frac{173}{14620} h^{11} f^{(10)}(\xi),

Анализ метода

Числовой пример

Рекомендации программисту

Автоматический выбор шага интегрирования с помощью апостериорной оценки погрешности методом Рунге

Величина погрешности численного интегрирования зависит как от шага сетки h, так и от гладкости подынтегральной функции f(x). В величину погрешности, кроме h входит также величина f^{(p)}(\xi), которая может сильно меняться на отрезке [a,b] и заранее неизвестна. Для уменьшения величины погрешности, можно измельчить сетку на заданном отрезке. Но при этом необходимо апостериорно оценивать погрешность. Такую оценку погрешности можно осуществить методом Рунге. Рассмотрим применение какой-то квадратурной формулы на частичном отрезке [x_{i-1},x_i]. Обозначим за I значение интеграла на всем отреpке, за I_i значение интеграла на i-ом частичном отрезке, за I_h приближенное значение интеграла на всем отрезке, полученное при помощи заданной квадратурной формулы и равномерном сетки с шагом h, а за I_{h,i} приближенное значение интеграла на i-ом частичном отрезке. Пусть данная квадратурная формула на данном частичном отрезке имеет порядок точности m, то есть

I_i-I_{h,i} \approx ch^m,

где с c - некоторая константа. Тогда

I_i-I_{ \frac{h}{2},i} \approx c( \frac{h}{2} )^m,

откуда получаем

I_i-I_{h,i} \approx 2^m(I_i-I_{ \frac{h}{2},i})
( 8 )
I_i-I_{ \frac{h}{2},i} \approx \frac {I_{ \frac{h}{2},i} - I_{h,i}} {2^m-1}.

Пусть используется составная квадратурная формула

I \approx I_h=\sum_{i=1}^N{I_{h,i}},

причем на всех частичных отрезках используются квадратурные формулы с одним и тем же порядком точности (или же,в частности, одна и та же формула). Проведем на каждом частичном отрезке [x_{i-1},x_i] вычисления дважды - один раз с шагом h_i, второй раз с шагом \frac{h_i}{2} и апостериорно оценим погрешность по правилу Рунге ( 8 ). Если для заданного \varepsilon >0 при i=1,2, \ldots N будут выполняться неравенства

( 9 )
|I_i-I_{ \frac{h}{2},i}| \approx \frac {|I_{ \frac{h}{2},i} - I_{h,i}|} {2^m-1} \le \frac {\varepsilon h_i} {b-a},

то получим

|I-I_{\frac{h}{2}}| \le \frac {\varepsilon}{b-a} \sum_{i=1}^N{h_i} \le \varepsilon,

то есть будет достигнута заданная точность \varepsilon.

Если же на каком-то из частичных отрезков оценка ( 9 ) не будет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида ( 9 ). Заметим, что для некоторой функции f(x) такое измельчение может продолжаться слишком долго. Поэтому в соответствующей программе следует предусмотреть ограничение сверху на число измельчений.

Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции f(x) и с мелким шагом - на участках быстрого изменения f(x). Это позволяет при заданной точности \epsilon уменьшить количество вычислений значений f(x) по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем, что для нахождения сумм I_{ \frac{h}{2},i} не надо пересчитывать значения f(x) во всех узлах, достаточно вычислять f(x) только в новых узлах.

Заключение

Список литературы

Личные инструменты