Методы парабол (Симпсона) и более высоких степеней (Ньютона - Котеса)
Материал из MachineLearning.
(9 промежуточных версий не показаны.) | |||
Строка 1: | Строка 1: | ||
== Введение == | == Введение == | ||
+ | |||
+ | В данной статье описывается метод вычисления собственного интеграла гладкой функции при помощи квадратурных формул. Формулы Ньютона-Котеса обладают следующими особенностями: | ||
+ | |||
+ | * Необходимым условием сходимости данного метода является существование и ограниченность производной функции (порядок производной зависит от выбранной формулы) | ||
+ | * Формулы Ньютона-Котеса обладают высоким порядком точности | ||
+ | * Формулы порядка <tex>n<9</tex> являются численно устойчивыми | ||
=== Постановка математической задачи === | === Постановка математической задачи === | ||
- | + | Задача численного интегрирования состоит в нахождении приближенного значения интеграла | |
{{ eqno | 1 }} | {{ eqno | 1 }} | ||
Строка 13: | Строка 19: | ||
где <tex>f(x_i)</tex> - значения функции <tex>f(x)</tex> в узлах <tex>x=x_i</tex> , где <tex>c_i</tex> - ''весовые множители'', зависящие только от узлов, но не зависящие от выбора <tex>f(x)</tex>. Формула {{eqref|2}} называется ''квадратурной формулой''. | где <tex>f(x_i)</tex> - значения функции <tex>f(x)</tex> в узлах <tex>x=x_i</tex> , где <tex>c_i</tex> - ''весовые множители'', зависящие только от узлов, но не зависящие от выбора <tex>f(x)</tex>. Формула {{eqref|2}} называется ''квадратурной формулой''. | ||
- | + | Задача численного интегрирования при помощи квадратур состоит в отыскании таких узлов <tex>\{x_i\}</tex> и таких весов <tex>\{c_i\}</tex>, чтобы ''погрешность квадратурной формулы'' | |
::<tex>D[f]=\sum_{i=0}^N{c_i f(x_i)} - \int_a^b{f(x)dx} = J_N[f] - J[f]</tex> | ::<tex>D[f]=\sum_{i=0}^N{c_i f(x_i)} - \int_a^b{f(x)dx} = J_N[f] - J[f]</tex> | ||
Строка 23: | Строка 29: | ||
::<tex>\int_a^b{f(x)dx}=\sum_{i=1}^N{\int_{x_{i-1}}^{x_i}{f(x)dx}}.</tex> | ::<tex>\int_a^b{f(x)dx}=\sum_{i=1}^N{\int_{x_{i-1}}^{x_i}{f(x)dx}}.</tex> | ||
- | Для построения формулы численного интегрирования на | + | Для построения формулы численного интегрирования на всём отрезке <tex>[a,b]</tex> достаточно построить квадратурную формулу для интеграла |
{{ eqno | 4 }} | {{ eqno | 4 }} | ||
Строка 32: | Строка 38: | ||
=== Построение квадратурных формул === | === Построение квадратурных формул === | ||
- | + | В силу вышеизложенного, вычисление приближенного значения интеграла производится при помощи квадратурной формулы | |
::<tex>\int_a^b{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.</tex> | ::<tex>\int_a^b{f(s)ds}\approx\sum_{k=0}^m{c_kf(s_k)}.</tex> | ||
Строка 60: | Строка 66: | ||
::<tex>c_0 s_0^m + c_1 s_1^m + \ldots + c_m s_m^m=\frac{1}{m+1}.</tex> | ::<tex>c_0 s_0^m + c_1 s_1^m + \ldots + c_m s_m^m=\frac{1}{m+1}.</tex> | ||
- | Эта система имеет единственное решение, так как | + | Эта система имеет единственное решение, так как её определителем является определитель Вандермонда, отличный от нуля если нет совпадающих узлов, <tex>s_0<s_q<\ldots<s_m</tex>. |
Так, полагая <tex>m=2,s_0=0,s_1=\frac{1}{2},s_2=1</tex>, имеем систему <tex>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}</tex>, решением которой являются веса '''формулы Симпсона''': <tex>c_0=c_2=\frac{1}{6},c_1=\frac{4}{6}</tex>. Таким образом, формула Симпсона является точной для полинома второй степени. Однако, в силу симметрии, она является точной и для всех полиномов третьей степени: | Так, полагая <tex>m=2,s_0=0,s_1=\frac{1}{2},s_2=1</tex>, имеем систему <tex>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}</tex>, решением которой являются веса '''формулы Симпсона''': <tex>c_0=c_2=\frac{1}{6},c_1=\frac{4}{6}</tex>. Таким образом, формула Симпсона является точной для полинома второй степени. Однако, в силу симметрии, она является точной и для всех полиномов третьей степени: | ||
Строка 85: | Строка 91: | ||
Формулы такого типа называют квадратурными '''формулами Котеса'''. | Формулы такого типа называют квадратурными '''формулами Котеса'''. | ||
+ | |||
+ | == Изложение метода == | ||
+ | |||
+ | === Применение квадратурных формул === | ||
+ | |||
+ | Вернемся к рассмотрению интеграла {{eqref|1}}. Как было показано выше, данный интеграл заменой сводится к интегралу на единичном отрезке <tex>[0,1]</tex>, а следовательно легко обобщить формулы для приближенного вычисления интеграла на единичном отрезке на произвольный. Применим аппарат квадратурных формул. Пусть задано равномерное разбиение отрезка с шагом <tex>h:</tex> <tex>\omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}</tex>, где обозначим <tex>f_i=f(a+ih)</tex>. Пусть также выбрана некоторая квадратурная формула Ньютона-Котеса (то есть выбрана степень полинома <tex>n</tex>, а следовательно каждый полином строится по <tex>n+1</tex> точке сетки). Мы также считаем, что данный набор из <tex>N+1</tex> точки можно разбить на поднаборы по <tex>m</tex> точек с совпадающими крайним, то есть | ||
+ | |||
+ | {{ eqno | 7 }} | ||
+ | ::<tex>N=n*k,</tex> где <tex>k \in \mathbb N.</tex> | ||
+ | |||
+ | Тогда, суммируя значения квадратур на каждом поднаборе, получим приближенное значение искомого интеграла. Если обозначить весовые множители <tex>c_i,i=0 \ldots n,</tex> то приближенное значение интеграла можно записать в виде двойной суммы | ||
+ | |||
+ | ::<tex>\int_a^b{f(x)dx \approx \sum_{j=0}^{k-1}{\sum_{i=0}^n{c_i f_{(i+j*n)}}}</tex> | ||
+ | |||
+ | Данный алгоритм естественным образом обобщается на случай, когда <tex>[a,b]=\bigcup_{i=0}^K{[a_i,b_i]}</tex>, где <tex>a_0=a, b_K=b, b_{i+1}=a_i;i=0\ldots K-1</tex>, а на каждом отрезке <tex>[a_i,b_i]</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> приближенное значение интеграла вычисляется при помощи квадратурных формул. | ||
=== Примеры квадратурных формул === | === Примеры квадратурных формул === | ||
- | + | Приведем примеры квадратурных формул Котеса на равномерной сетке с шагом <tex>h:</tex> <tex>\omega_h=\{x_i=a+ih, i=0,1,\ldots,N;hN=b-a}</tex>, где обозначим <tex>f_i=f(a+ih),</tex> <tex>\xi \in [a,b]</tex>: | |
+ | |||
+ | *для 2 точек([[Методы прямоугольников и трапеций|метод трапеций]]), <tex>m=1:</tex> | ||
+ | |||
+ | ::<tex>J_1[f]=\frac{1}{2}h(f_0+f_1),</tex> | ||
+ | ::<tex>D[f]=\frac{1}{12} h^3 f^{(2)}(\xi),</tex> | ||
+ | |||
*для 3 точек(метод Симпсона), <tex>m=2:</tex> | *для 3 точек(метод Симпсона), <tex>m=2:</tex> | ||
::<tex>J_2[f]=\frac{1}{3}h(f_0+4f_1+f_2),</tex> | ::<tex>J_2[f]=\frac{1}{3}h(f_0+4f_1+f_2),</tex> | ||
- | ::<tex> | + | ::<tex>D[f]=\frac{1}{90} h^5 f^{(4)}(\xi),</tex> |
*для 4 точек, <tex>m=3:</tex> | *для 4 точек, <tex>m=3:</tex> | ||
::<tex>J_3[f]=\frac{3}{8}h(f_0+3f_1+3f_2+f_3),</tex> | ::<tex>J_3[f]=\frac{3}{8}h(f_0+3f_1+3f_2+f_3),</tex> | ||
- | ::<tex> | + | ::<tex>D[f]=\frac{3}{80} h^5 f^{(4)}(\xi),</tex> |
*для 5 точек, <tex>m=4:</tex> | *для 5 точек, <tex>m=4:</tex> | ||
::<tex>J_4[f]=\frac{2}{45}h(7f_0+32f_1+12f_2+32f_3+7f_4),</tex> | ::<tex>J_4[f]=\frac{2}{45}h(7f_0+32f_1+12f_2+32f_3+7f_4),</tex> | ||
- | ::<tex> | + | ::<tex>D[f]=\frac{8}{945} h^7 f^{(6)}(\xi),</tex> |
*для 6 точек, <tex>m=5:</tex> | *для 6 точек, <tex>m=5:</tex> | ||
::<tex>J_5[f]=\frac{5}{288}h(19f_0+75f_1+50f_2+50f_3+75f_4+19f_5),</tex> | ::<tex>J_5[f]=\frac{5}{288}h(19f_0+75f_1+50f_2+50f_3+75f_4+19f_5),</tex> | ||
+ | ::<tex>D[f]=\frac{275}{12096} h^7 f^{(6)}(\xi),</tex> | ||
*для 7 точек, <tex>m=6:</tex> | *для 7 точек, <tex>m=6:</tex> | ||
::<tex>J_6[f]=\frac{1}{140}h(41f_0+216f_1+27f_2+272f_3+27f_4+216f_5+41f_6),</tex> | ::<tex>J_6[f]=\frac{1}{140}h(41f_0+216f_1+27f_2+272f_3+27f_4+216f_5+41f_6),</tex> | ||
+ | ::<tex>D[f]=\frac{9}{1400} h^9 f^{(8)}(\xi),</tex> | ||
*для 8 точек, <tex>m=7:</tex> | *для 8 точек, <tex>m=7:</tex> | ||
::<tex>J_7[f]=\frac{7}{17280}h(751f_0+3577f_1+1323f_2+2989f_3+2989f_4+1323f_5+3577f_6+751f_7),</tex> | ::<tex>J_7[f]=\frac{7}{17280}h(751f_0+3577f_1+1323f_2+2989f_3+2989f_4+1323f_5+3577f_6+751f_7),</tex> | ||
+ | ::<tex>D[f]=\frac{8183}{518400} h^9 f^{(8)}(\xi),</tex> | ||
*для 9 точек, <tex>m=8:</tex> | *для 9 точек, <tex>m=8:</tex> | ||
::<tex>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),</tex> | ::<tex>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),</tex> | ||
+ | ::<tex>D[f]=\frac{2368}{467775} h^{11} f^{(10)}(\xi),</tex> | ||
*для 10 точек, <tex>m=9:</tex> | *для 10 точек, <tex>m=9:</tex> | ||
::<tex>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).</tex> | ::<tex>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).</tex> | ||
- | + | ::<tex>D[f]=\frac{173}{14620} h^{11} f^{(10)}(\xi),</tex> | |
- | + | ||
== Анализ метода == | == Анализ метода == | ||
- | == Числовой | + | === Погрешность квадратурной формулы === |
+ | |||
+ | Пусть функция <tex>f</tex> имеет <tex>(n+1)</tex>-ю непрерывную производную на отрезке <tex>0 \le s \le 1</tex>, то есть <tex>f(s) \in C^{(n+1)}[0,1]</tex>, <tex>s_k</tex> - точки, в которых задано значение функции. Пусть используется квадратурная формула <tex>n</tex>-го порядка. Введем функцию | ||
+ | |||
+ | ::<tex>K_n(\xi)= \begin{cases} \frac {\xi ^n} {n!} ,& \xi \ge 0;\\ 0 ,& \xi <0. \end{cases}</tex> | ||
+ | |||
+ | Тогда формула для погрешности имеет следующий вид | ||
+ | |||
+ | ::<tex>D[f]=\int_0^1{F_{n+1}(t)f^{(n+1)}(t)dt},</tex> | ||
+ | |||
+ | где | ||
+ | |||
+ | ::<tex>F_{n+1}(t)=\sum_{k=0}^m {c_k K_n(s_k-t)} - \frac {(1-t)^{n+1}} {(n+1)!}</tex> | ||
+ | |||
+ | Отсюда следует оценка для погрешности | ||
+ | |||
+ | {{ eqno | 8 }} | ||
+ | ::<tex>|D[f]| \le M_{n+1} z_{n+1}</tex> | ||
+ | |||
+ | при <tex>|f^{(n+1)}(t)| \le M_{n+1}</tex>, где <tex>M_{n+1}>0</tex> - постоянная, и при | ||
+ | |||
+ | ::<tex>z_{n+1}=\int_0^1{|F_{n+1}(t)|dt}.</tex> | ||
+ | |||
+ | Если <tex>F_{n+1}(t)</tex> не меняет знака на отрезке <tex>0 \le s \le 1</tex>, то в силу теоремы о среднем имеем | ||
+ | |||
+ | {{ eqno | 9 }} | ||
+ | ::<tex>D[f]=f^{(n+1)}(\xi) \int_0^1 {F_{n+1}(t)dt}, \xi \in [0,1].</tex> | ||
+ | |||
+ | === Численная устойчивость квадратурных формул === | ||
+ | |||
+ | Отметим, что формулы Ньютона-Котеса с <tex>n \ge 10</tex> редко используются из-за их численной неустойчивости, приводящей к резкому возрастанию вычислительной погрешности. Причиной такой неустойчивости является то, что коэффициенты формул Ньютона-Котеса при больших <tex>n</tex> имеют различные знаки, а именно при <tex>n \ge 10</tex> существуют как положительные, так и отрицательные коэффициенты. | ||
+ | |||
+ | Рассмотрим квадратурную сумму | ||
+ | |||
+ | {{ eqno | 10 }} | ||
+ | ::<tex>I_n= \sum_{k=0}^n {c_k f(x_k)}.</tex> | ||
+ | |||
+ | Предположим, что значения функции <tex>f(x)</tex>, заданной на отрезке <tex>[a,b]</tex>, вычисляются с некоторой погрешностью, то есть вместо точного значения получаем приближенное значение <tex>\bar f(x_k)=f(x_k)+ \delta _k</tex>. Тогда вместо <tex>I_n</tex> получим сумму | ||
+ | |||
+ | ::<tex>\bar I_n= \sum_{k=0}^n {c_k(f(x_k)+\delta _k)}=I_n + \delta I_n,</tex> | ||
+ | |||
+ | где | ||
+ | |||
+ | {{ eqno | 11 }} | ||
+ | ::<tex>\delta I_n=\sum_{k=0}^n {c_k \delta _k}.</tex> | ||
+ | |||
+ | Поскольку квадратурная формула точна для <tex>f(x) \equiv 1</tex>, имеем | ||
+ | |||
+ | {{ eqno | 12 }} | ||
+ | ::<tex>\sum_{k=0}^n{c_k}=\int_a^b{dx}=b-a</tex> | ||
+ | |||
+ | и не зависит от <tex>n</tex>. | ||
+ | |||
+ | Предположим теперь, что все коэффициенты <tex>c_k</tex> неотрицательны. Тогда из {{ eqref | 11 }} и {{ eqref | 12 }} получим оценку | ||
+ | |||
+ | {{ eqno | 13 }} | ||
+ | ::<tex>|\delta I_n| \le \sum_{k=0}^n {|c_k||\delta _k|}=\sum_{k=0}^n{c_k|\delta _k|} \le (\mathrm{}\max_{0 \le k \le n} |\delta _k|)(b-a),</tex> | ||
+ | |||
+ | которая означает, что при больших <tex>n</tex> погрешность в вычислении квадратурной суммы {{ eqref | 10 }} имеет тот же порядок, что и погрешность в вычислении функции. В этом случае говорят, что сумма {{ eqref | 10 }} вычисляется устойчиво. | ||
+ | |||
+ | Если коэффициенты <tex>c_k</tex> имеют различные знаки, то может оказаться, что сумма <tex>\sum_{k=0}^n{|c_k|}</tex> не является равномерно ограниченной по <tex>n</tex> и, следовательно погрешность в вычислении <tex>I_n</tex> неограниченно возрастает с ростом <tex>n</tex>. В этом случае вычисления по формуле {{ eqref | 10 }} будут неустойчивы и пользоваться такой формулой при больших <tex>n</tex> нельзя. | ||
+ | |||
+ | == Числовой эксперимент == | ||
+ | |||
+ | Приведем примеры вычисления интегралов с использованием формул Ньютона-Котеса. При реализации использовался язык C++, ниже приведен код функции, возвращающей приближенное значение интеграла. | ||
+ | |||
+ | === Исходный код функции === | ||
+ | |||
+ | На вход функция принимает 4 параметра | ||
+ | |||
+ | double a - левый конец исследуемого отрезка | ||
+ | |||
+ | double b - правый конец исследуемого отрезка | ||
+ | |||
+ | int Degree - степень используемого полинома | ||
+ | |||
+ | int Ndivisions - количество отрезков, на которые разбивается исходный. Совпадает со значением <tex>k</tex> в формуле {{ eqref | 7 }} | ||
+ | |||
+ | f - интегрируемая функция | ||
+ | |||
+ | <pre> | ||
+ | double f(double x) | ||
+ | { | ||
+ | return 1/x; | ||
+ | } | ||
+ | |||
+ | double NewtonCotes(double a, double b, int Degree, int Ndivisions) | ||
+ | { | ||
+ | int koef[10][10]={ 1,0,0,0,0,0,0,0,0,0, | ||
+ | 1,1,0,0,0,0,0,0,0,0, | ||
+ | 1,4,1,0,0,0,0,0,0,0, | ||
+ | 1,3,3,1,0,0,0,0,0,0, | ||
+ | 7,32,12,32,7,0,0,0,0,0, | ||
+ | 19,75,50,50,75,19,0,0,0,0, | ||
+ | 41,216,27,272,27,216,41,0,0,0, | ||
+ | 751,3577,1323,2989,2989,1323,3577,751,0,0, | ||
+ | 989,5888,-928,10496,-4540,10496,-928,5888,989,0, | ||
+ | 2857,15741,1080,19344,5778,5778,19344,1080,15741,2857 | ||
+ | }; | ||
+ | double mltp[10]={1,1.0/2,1.0/3,3.0/8,2.0/45,5.0/288,1.0/140,7.0/17280,4.0/14175,9.0/89600}; | ||
+ | |||
+ | if ((Degree<0) || (Degree>9))throw "Wrong degree"; | ||
+ | if (a>=b) throw "Wrong segment"; | ||
+ | if (Ndivisions<1) Ndivisions = 1; | ||
+ | |||
+ | double Sum,PartSum; | ||
+ | double h=(b-a)/(Degree*Ndivisions); | ||
+ | |||
+ | Sum=0; | ||
+ | for (int j=0; j<Ndivisions; j++) | ||
+ | { | ||
+ | PartSum=0; | ||
+ | for (int i=0; i<=Degree; i++) | ||
+ | PartSum+= koef[Degree][i]*f(a + (i+j*Degree) * h ); | ||
+ | Sum+= mltp[Degree]*PartSum*h; | ||
+ | } | ||
+ | |||
+ | return Sum; | ||
+ | } | ||
+ | </pre> | ||
+ | |||
+ | === Пример === | ||
+ | |||
+ | Вычислим по формулам Котеса степеней от 1 до 9 значение интеграла | ||
+ | |||
+ | ::<tex>I=\int_1^2{ \frac{dx}{x} } = ln(2)</tex> | ||
+ | |||
+ | Вычисления будем производить не разбивая отрезок на частичные, то есть используя <tex>k+1</tex> точку, где <tex>k</tex>-степень полинома. Результаты приведены ниже в таблице. Округления проводились с точностью 6 знаков после запятой. | ||
+ | |||
+ | {| class="wikitable" | ||
+ | |- | ||
+ | ! Степень полинома | ||
+ | ! Количество узлов | ||
+ | ! Приближенное значение интеграла | ||
+ | ! Точное значение интеграла | ||
+ | ! Погрешность результата | ||
+ | ! Погрешность формулы | ||
+ | |- | ||
+ | | 1 | ||
+ | | 2 | ||
+ | | 0.75 | ||
+ | | 0.693147 | ||
+ | | 0.056853 | ||
+ | | 0.166667 | ||
+ | |- | ||
+ | | 2 | ||
+ | | 3 | ||
+ | | 0.694444 | ||
+ | | 0.693147 | ||
+ | | 0.001297 | ||
+ | | 0.008333 | ||
+ | |- | ||
+ | | 3 | ||
+ | | 4 | ||
+ | | 0.69375 | ||
+ | | 0.693147 | ||
+ | | 0.000603 | ||
+ | | 0.003704 | ||
+ | |- | ||
+ | | 4 | ||
+ | | 5 | ||
+ | | 0.693175 | ||
+ | | 0.693147 | ||
+ | | 0.000028 | ||
+ | | 0.000372 | ||
+ | |- | ||
+ | | 5 | ||
+ | | 6 | ||
+ | | 0.693163 | ||
+ | | 0.693147 | ||
+ | | 0.000016 | ||
+ | | 0.000210 | ||
+ | |- | ||
+ | | 6 | ||
+ | | 7 | ||
+ | | 0.693148 | ||
+ | | 0.693147 | ||
+ | | 0.000001 | ||
+ | | 0.000026 | ||
+ | |- | ||
+ | | 7 | ||
+ | | 8 | ||
+ | | 0.693148 | ||
+ | | 0.693147 | ||
+ | | 0.000001 | ||
+ | | 0.000016 | ||
+ | |- | ||
+ | | 8 | ||
+ | | 9 | ||
+ | | 0.693147 | ||
+ | | 0.693147 | ||
+ | | 0 | ||
+ | | 0.000002 | ||
+ | |- | ||
+ | | 9 | ||
+ | | 10 | ||
+ | | 0.693147 | ||
+ | | 0.693147 | ||
+ | | 0 | ||
+ | | 0 | ||
+ | |} | ||
== Рекомендации программисту == | == Рекомендации программисту == | ||
+ | |||
+ | === Автоматический выбор шага интегрирования с помощью апостериорной оценки погрешности методом Рунге === | ||
+ | |||
+ | Величина погрешности численного интегрирования зависит как от шага сетки <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 | 14 }} | ||
+ | ::<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 | 14 }}. Если для заданного <tex>\varepsilon >0</tex> при <tex>i=1,2, \ldots N</tex> будут выполняться неравенства | ||
+ | |||
+ | {{eqno | 15 }} | ||
+ | ::<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 | 15 }} не будет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида {{ eqref | 15 }}. Заметим, что для некоторой функции <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> только в новых узлах. | ||
== Заключение == | == Заключение == | ||
+ | |||
+ | Формулы Симпсона и Ньютона-Котеса являются хорошим аппаратом для вычисления определенного интеграла достаточное число раз непрерывно дифференцируемой функции. На примере формул <tex>n</tex>-го порядка мы видим, что при ограниченной производной <tex>(n+1)</tex>-го порядка точность формул есть <tex>O(h^p)</tex>, где <tex>h</tex> - шаг сетки, а <tex>p>n</tex>. То есть при выполнении условий применимости данного метода, формулы дают высокий порядок сходимости. Однако при больших <tex>n</tex>, в частности уже при <tex>n \ge 10</tex>, вычисление приближенного значения интеграла становится численно неустойчиво, что делает их неприменимыми. | ||
== Список литературы == | == Список литературы == | ||
Строка 134: | Строка 406: | ||
* http://mathworld.wolfram.com/Newton-CotesFormulas.html | * http://mathworld.wolfram.com/Newton-CotesFormulas.html | ||
- | + | == См. также == | |
+ | *[[Практикум ММП ВМК, 4й курс, осень 2008]] | ||
+ | *[[Media:Newton_Cotes_Example1.zip| Код программы ]] | ||
+ | |||
[[Категория:Численное интегрирование]] | [[Категория:Численное интегрирование]] | ||
+ | [[Категория:Учебные задачи]] |
Текущая версия
Содержание |
Введение
В данной статье описывается метод вычисления собственного интеграла гладкой функции при помощи квадратурных формул. Формулы Ньютона-Котеса обладают следующими особенностями:
- Необходимым условием сходимости данного метода является существование и ограниченность производной функции (порядок производной зависит от выбранной формулы)
- Формулы Ньютона-Котеса обладают высоким порядком точности
- Формулы порядка являются численно устойчивыми
Постановка математической задачи
Задача численного интегрирования состоит в нахождении приближенного значения интеграла
где - заданная и интегрируемая на отрезке функция. На отрезке вводится сетка и в качестве приближенного значения интеграла рассматривается число
где - значения функции в узлах , где - весовые множители, зависящие только от узлов, но не зависящие от выбора . Формула (2) называется квадратурной формулой. Задача численного интегрирования при помощи квадратур состоит в отыскании таких узлов и таких весов , чтобы погрешность квадратурной формулы
была минимальной по модулю для функции из заданного класса (величина зависит от гладкости ). Погрешность зависит как от расположения узлов, так и от выбора весовых коэффициентов. Введем на равномерную сетку с шагом , т.е. множество точек , и представим интеграл (1) в виде суммы интегралов по частичным отрезкам:
Для построения формулы численного интегрирования на всём отрезке достаточно построить квадратурную формулу для интеграла
на частичном отрезке и воспользоваться свойством (3).
Построение квадратурных формул
В силу вышеизложенного, вычисление приближенного значения интеграла производится при помощи квадратурной формулы
Данную формулу при помощи замены можно привести к стандартному виду
В общем случае узлы и веса неизвестны и подлежат определению.
Рассмотрим случай, когда узлы заданы и требуется найти веса квадратурной формулы . Будем пользоваться требованием: формула (5) должна быть точной для любого полинома степени . Для того чтобы полином степени удовлетворял данному требованию, достаточно потребовать, чтобы квадратурная формула была точной для любого одночлена степени . Учитывая, что , получаем уравнение
Эта система имеет единственное решение, так как её определителем является определитель Вандермонда, отличный от нуля если нет совпадающих узлов, .
Так, полагая , имеем систему , решением которой являются веса формулы Симпсона: . Таким образом, формула Симпсона является точной для полинома второй степени. Однако, в силу симметрии, она является точной и для всех полиномов третьей степени:
так как она точна для :
Формулы треугольника и трапеции точны для линейной функции,т.е. для полинома первой степени, в чем легко убедиться непосредственно. В общем случае в качестве можно выбрать интерполяционный полином Лагранжа
где - интерполяционный коэффициент Лагранжа. Из равенства
видно, что формула (5) верна для полинома степени , если весовые коэффициенты определяются по формуле
Формулы такого типа называют квадратурными формулами Котеса.
Изложение метода
Применение квадратурных формул
Вернемся к рассмотрению интеграла (1). Как было показано выше, данный интеграл заменой сводится к интегралу на единичном отрезке , а следовательно легко обобщить формулы для приближенного вычисления интеграла на единичном отрезке на произвольный. Применим аппарат квадратурных формул. Пусть задано равномерное разбиение отрезка с шагом , где обозначим . Пусть также выбрана некоторая квадратурная формула Ньютона-Котеса (то есть выбрана степень полинома , а следовательно каждый полином строится по точке сетки). Мы также считаем, что данный набор из точки можно разбить на поднаборы по точек с совпадающими крайним, то есть
- где
Тогда, суммируя значения квадратур на каждом поднаборе, получим приближенное значение искомого интеграла. Если обозначить весовые множители то приближенное значение интеграла можно записать в виде двойной суммы
Данный алгоритм естественным образом обобщается на случай, когда , где , а на каждом отрезке задана равномерная сетка. Тогда искомый интеграл равен
а на каждом из частичных отрезков приближенное значение интеграла вычисляется при помощи квадратурных формул.
Примеры квадратурных формул
Приведем примеры квадратурных формул Котеса на равномерной сетке с шагом , где обозначим :
- для 2 точек(метод трапеций),
- для 3 точек(метод Симпсона),
- для 4 точек,
- для 5 точек,
- для 6 точек,
- для 7 точек,
- для 8 точек,
- для 9 точек,
- для 10 точек,
Анализ метода
Погрешность квадратурной формулы
Пусть функция имеет -ю непрерывную производную на отрезке , то есть , - точки, в которых задано значение функции. Пусть используется квадратурная формула -го порядка. Введем функцию
Тогда формула для погрешности имеет следующий вид
где
Отсюда следует оценка для погрешности
при , где - постоянная, и при
Если не меняет знака на отрезке , то в силу теоремы о среднем имеем
Численная устойчивость квадратурных формул
Отметим, что формулы Ньютона-Котеса с редко используются из-за их численной неустойчивости, приводящей к резкому возрастанию вычислительной погрешности. Причиной такой неустойчивости является то, что коэффициенты формул Ньютона-Котеса при больших имеют различные знаки, а именно при существуют как положительные, так и отрицательные коэффициенты.
Рассмотрим квадратурную сумму
Предположим, что значения функции , заданной на отрезке , вычисляются с некоторой погрешностью, то есть вместо точного значения получаем приближенное значение . Тогда вместо получим сумму
где
Поскольку квадратурная формула точна для , имеем
и не зависит от .
Предположим теперь, что все коэффициенты неотрицательны. Тогда из ( 11 ) и ( 12 ) получим оценку
которая означает, что при больших погрешность в вычислении квадратурной суммы ( 10 ) имеет тот же порядок, что и погрешность в вычислении функции. В этом случае говорят, что сумма ( 10 ) вычисляется устойчиво.
Если коэффициенты имеют различные знаки, то может оказаться, что сумма не является равномерно ограниченной по и, следовательно погрешность в вычислении неограниченно возрастает с ростом . В этом случае вычисления по формуле ( 10 ) будут неустойчивы и пользоваться такой формулой при больших нельзя.
Числовой эксперимент
Приведем примеры вычисления интегралов с использованием формул Ньютона-Котеса. При реализации использовался язык C++, ниже приведен код функции, возвращающей приближенное значение интеграла.
Исходный код функции
На вход функция принимает 4 параметра
double a - левый конец исследуемого отрезка
double b - правый конец исследуемого отрезка
int Degree - степень используемого полинома
int Ndivisions - количество отрезков, на которые разбивается исходный. Совпадает со значением в формуле ( 7 )
f - интегрируемая функция
double f(double x) { return 1/x; } double NewtonCotes(double a, double b, int Degree, int Ndivisions) { int koef[10][10]={ 1,0,0,0,0,0,0,0,0,0, 1,1,0,0,0,0,0,0,0,0, 1,4,1,0,0,0,0,0,0,0, 1,3,3,1,0,0,0,0,0,0, 7,32,12,32,7,0,0,0,0,0, 19,75,50,50,75,19,0,0,0,0, 41,216,27,272,27,216,41,0,0,0, 751,3577,1323,2989,2989,1323,3577,751,0,0, 989,5888,-928,10496,-4540,10496,-928,5888,989,0, 2857,15741,1080,19344,5778,5778,19344,1080,15741,2857 }; double mltp[10]={1,1.0/2,1.0/3,3.0/8,2.0/45,5.0/288,1.0/140,7.0/17280,4.0/14175,9.0/89600}; if ((Degree<0) || (Degree>9))throw "Wrong degree"; if (a>=b) throw "Wrong segment"; if (Ndivisions<1) Ndivisions = 1; double Sum,PartSum; double h=(b-a)/(Degree*Ndivisions); Sum=0; for (int j=0; j<Ndivisions; j++) { PartSum=0; for (int i=0; i<=Degree; i++) PartSum+= koef[Degree][i]*f(a + (i+j*Degree) * h ); Sum+= mltp[Degree]*PartSum*h; } return Sum; }
Пример
Вычислим по формулам Котеса степеней от 1 до 9 значение интеграла
Вычисления будем производить не разбивая отрезок на частичные, то есть используя точку, где -степень полинома. Результаты приведены ниже в таблице. Округления проводились с точностью 6 знаков после запятой.
Степень полинома | Количество узлов | Приближенное значение интеграла | Точное значение интеграла | Погрешность результата | Погрешность формулы |
---|---|---|---|---|---|
1 | 2 | 0.75 | 0.693147 | 0.056853 | 0.166667 |
2 | 3 | 0.694444 | 0.693147 | 0.001297 | 0.008333 |
3 | 4 | 0.69375 | 0.693147 | 0.000603 | 0.003704 |
4 | 5 | 0.693175 | 0.693147 | 0.000028 | 0.000372 |
5 | 6 | 0.693163 | 0.693147 | 0.000016 | 0.000210 |
6 | 7 | 0.693148 | 0.693147 | 0.000001 | 0.000026 |
7 | 8 | 0.693148 | 0.693147 | 0.000001 | 0.000016 |
8 | 9 | 0.693147 | 0.693147 | 0 | 0.000002 |
9 | 10 | 0.693147 | 0.693147 | 0 | 0 |
Рекомендации программисту
Автоматический выбор шага интегрирования с помощью апостериорной оценки погрешности методом Рунге
Величина погрешности численного интегрирования зависит как от шага сетки , так и от гладкости подынтегральной функции . В величину погрешности, кроме входит также величина , которая может сильно меняться на отрезке и заранее неизвестна. Для уменьшения величины погрешности, можно измельчить сетку на заданном отрезке. Но при этом необходимо апостериорно оценивать погрешность. Такую оценку погрешности можно осуществить методом Рунге. Рассмотрим применение какой-то квадратурной формулы на частичном отрезке . Обозначим за значение интеграла на всем отреpке, за значение интеграла на -ом частичном отрезке, за приближенное значение интеграла на всем отрезке, полученное при помощи заданной квадратурной формулы и равномерном сетки с шагом , а за приближенное значение интеграла на -ом частичном отрезке. Пусть данная квадратурная формула на данном частичном отрезке имеет порядок точности , то есть
где с - некоторая константа. Тогда
откуда получаем
Пусть используется составная квадратурная формула
причем на всех частичных отрезках используются квадратурные формулы с одним и тем же порядком точности (или же,в частности, одна и та же формула). Проведем на каждом частичном отрезке вычисления дважды - один раз с шагом , второй раз с шагом и апостериорно оценим погрешность по правилу Рунге ( 14 ). Если для заданного при будут выполняться неравенства
то получим
то есть будет достигнута заданная точность .
Если же на каком-то из частичных отрезков оценка ( 15 ) не будет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида ( 15 ). Заметим, что для некоторой функции такое измельчение может продолжаться слишком долго. Поэтому в соответствующей программе следует предусмотреть ограничение сверху на число измельчений.
Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции и с мелким шагом - на участках быстрого изменения . Это позволяет при заданной точности уменьшить количество вычислений значений по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем, что для нахождения сумм не надо пересчитывать значения во всех узлах, достаточно вычислять только в новых узлах.
Заключение
Формулы Симпсона и Ньютона-Котеса являются хорошим аппаратом для вычисления определенного интеграла достаточное число раз непрерывно дифференцируемой функции. На примере формул -го порядка мы видим, что при ограниченной производной -го порядка точность формул есть , где - шаг сетки, а . То есть при выполнении условий применимости данного метода, формулы дают высокий порядок сходимости. Однако при больших , в частности уже при , вычисление приближенного значения интеграла становится численно неустойчиво, что делает их неприменимыми.
Список литературы
- А.А.Самарский, А.В.Гулин. Численные методы М.: Наука, 1989.
- А.А.Самарский. Введение в численные методы М.: Наука, 1982.
- http://mathworld.wolfram.com/Newton-CotesFormulas.html