\documentclass[12pt,a4paper,oneside,reqno,draft]{amsbook} \usepackage{amsmath} \usepackage{amsthm} \usepackage[T2A]{fontenc} \usepackage[cp1251]{inputenc} \usepackage[english,russian]{babel} \usepackage{amsfonts} \usepackage{latexsym} \tolerance 900 \newtheorem{lemma}{Лемма} \newtheorem{theorem}{Теорема}[section] \newtheorem{example}{Пример}[section] \newtheorem{proposition}{Предложение} \newtheorem{corollary}{Следствие} \makeatletter \@addtoreset{equation}{section} \makeatother \renewcommand{\theequation}{\arabic{section}.\arabic{equation}} \renewcommand{\thetheorem}{\arabic{section}.\arabic{theorem}} \begin{document} \bigskip \begin{center}\bf Государственный комитет Российской Федерации\\ по высшему образованию \vskip20pt \large Московский государственный авиационный\\ технологический университет им.\ К.~Э.~Циолковского \bigskip \bigskip Кафедра ``Высшая математика'' \vskip100pt \vskip70pt \Large КВАДРАТУРНЫЕ ФОРМУЛЫ \vskip50pt \rm Методические указания по курсу ``Численные методы'' \vskip70pt \begin{flushright} Составитель: Осипенко К.Ю. \end{flushright} \vskip 90pt Москва 1995 г. \end{center} \thispagestyle{empty} \pagestyle{plain} \section*{\bf Введение} Если $f(x)$ --- непрерывная на отрезке $[a,b]$ функция и $F(x)$ --- ее первообразная, то по формуле Ньютона--Лейбница $$\int_a^bf(x)\,dx=F(a)-F(b).$$ Однако часто первообразная не может быть выражена через элементарные функции или является слишком сложной. Например, при вычислении интегралов $$\int_0^1e^{-x^2}\,dx,\qquad\int_0^1\frac{dx}{(1+x^2)^{10}}$$ в первом случае мы сталкиваемся с тем, что первообразная для $e^{-x^2}$ не выражается через элементарные функции, а во втором --- с тем, что первообразная для функции $(1+x^2)^{-10}$ является слишком громоздким выражением. Кроме того, на практике подынтегральная функция $f(x)$ обычно бывает задана в дискретном числе точек. В этом случае первообразная $F(x)$ вообще не может быть найдена точно. Тем самым возникает задача приближенного вычисления определенного интеграла от функции по информации о значениях этой функции в некоторой системе точек. Такого рода формулы называются {\it квадратурными формулами}. В данном пособии рассматриваются основные методы построения простейших квадратурных формул и оценки их погрешностей. \section{Формула прямоугольников} Рассмотрим задачу приближенного вычисления интеграла $$\int_{-h/2}^{h/2}f(x)\,dx$$ по значению $f(0)=f_0$. Естественно считать, что функция всюду приближенно равна $f_0$, и заменить вычисление интеграла от исходной функции на вычисление интеграла от постоянной $f_0$. Таким образом, мы приходим к простейшей квадратурной формуле $$\int_{-h/2}^{h/2}f(x)\,dx\approx hf_0,$$ называемой {\it формулой прямоугольников}. Геометрический смысл формулы прямоугольников заключается в том, что площадь под графиком функции $y=f(x)$ заменяется на площадь прямоугольника с высотой, равной $f(0)$. Оценим погрешность формулы прямоугольников в предположении, что у функции $f( x)$ существует непрерывная вторая производная. По формуле Тейлора имеем $$f(x)=f_0+f_0'x+\frac{f''(\xi)}{2!}x^2,$$ где $f_0'=f'(0)$. Отсюда $$\int_{-h/2}^{h/2}f(x)\,dx=hf_0+\int_{-h/2}^{h/2}\frac{f''(\xi)}{2!}x^2\,dx. $$ Положим $$M_2=\max_{x\in[-h/2,h/2]}|f''(x)|.$$ Тогда \begin{equation}\label{1.1} \left|\int_{-h/2}^{h/2}f(x)\,dx-hf_0\right|\le\frac{M_2}2\int_{-h/2}^{h/2}x ^2\,dx=\frac{M_2}{24}h^3. \end{equation} \section{Усложненная формула прямоугольников} Пусть имеется отрезок $[a,b]$ и требуется приближенно вычислить $$\int_a^bf(x)\,dx.$$ Разобьем отрезок $[a,b]$ на $N$ равных частей точками $$x_i=a+ih,\quad i=0,1,\ldots,N,\quad h=\frac{b-a}N.$$ На каждом из отрезков $[x_i,x_{i+1}]$, $i=0,\ldots,N-1$, вычислим значение функции в средней точке $x_{i+1/2}=a+(i+1/2)h$ и применим формулу прямоугольников \begin{equation}\label{2.1} \int_{x_i}^{x_{i+1}}f(x)\,dx\approx hf_{i+1/2}, \end{equation} где $f_{i+1/2}=f(x_{i+1/2})$. Поскольку $$\int_a^bf(x)\,dx=\sum_{i=0}^{N-1}\int_{x_i}^{x_{i+1}}f(x)\,dx,$$ то, сложив приближенные равенства \eqref{2.1}, получим \begin{equation}\label{2.2} \int_a^bf(x)\,dx\approx h(f_{1/2}+f_{3/2}+\ldots+f_{N-1/2}). \end{equation} Формула \eqref{2.2} называется {\it усложненной формулой прямоугольников\/} (часто именно эту формулу называют {\it формулой прямоугольников\/}). Для оценки погрешности формулы \eqref{2.2} предположим, что функция $f(x)$ имеет непрерывную вторую производную на отрезке $[a,b]$ и $$M_2=\max_{x\in[a,b]}|f''(x)|.$$ Учитывая оценку \eqref{1.1}, имеем $$\biggl|\int_a^bf(x)\,dx-h\sum_{i=0}^{N-1}f_{i+1/2}\biggr|\le\sum_{i=0}^{N-1 }\frac{M_2}{24}h^3=\frac{M_2}{24}(b-a)h^2.$$ В силу равенства $Nh=b-a$ можно переписать эту оценку следующим образом $$\biggl|\int_a^bf(x)\,dx-h\sum_{i=0}^{N-1}f_{i+1/2}\biggr|\le\frac{M_2(b-a) ^3}{24N^2}.$$ \section{Использование интерполяционного многочлена Лагранжа для построения квадратурных формул} Одним из общих приемов построения квадратурных формул является замена функции, заданной на отрезке $[a,b]$, некоторой более простой и в то же время близкой к исходной функцией. Например, если $f(x)$ известна в некоторых точках отрезка $[a,b]$ \ $x_0,x_1,\ldots,x_n$, то можно заменить ее на интерполяционный многочлен Лагранжа $L_n(x)$ и положить $$\int_a^bf(x)\,dx\approx\int_a^bL_n(x)\,dx.$$ Для интерполяционного многочлена Лагранжа имеет место равенство (см.\ [7, стр.\ 9]) $$L_n(x)=\sum_{i=0}^nl_{ni}(x)f_i,$$ где $$l_{ni}(x)=\frac{(x-x_0)\ldots(x-x_{i-1})(x-x_{i+1})\ldots(x-x_n)}{(x_i-x_0) \ldots(x_i-x_{i-1})(x_i-x_{i+1})\ldots(x_i-x_n)},\quad f_i=f(x_i).$$ Таким образом, получаем квадратурную формулу \begin{equation}\label{3.1} \int_a^bf(x)\,dx\approx\sum_{i=0}^nA_if_i, \end{equation} в которой \begin{equation}\label{3.2} A_i=\int_a^bl_{ni}(x)\,dx. \end{equation} {\sc Замечание.} Если $f(x)$ --- многочлен степени $k\le n$, то в силу единственности интерполяционного многочлена Лагранжа $f(x)\equiv L_n(x)$. Тем самым квадратурная формула \eqref{3.1} точна на многочленах степени $k\le n$. В частности, при всех $k=0,1,\ldots,n$ \begin{equation}\label{3.3} \int_a^bx^k\,dx=\sum_{i=0}^nA_ix_i^k. \end{equation} \begin{example}\rm Построить квадратурную формулу вида $$\int_{-1}^1f(x)\,dx\approx A_0f(-1/2)+A_1f(0)+A_2f(1/4).$$ \end{example} {\sc Решение.} Можно было бы вычислять $A_0$, $A_1$ и $A_2$ по формулам \eqref{3.2}. Например, $$A_0=\int_{-1}^1\frac{x(x-1/4)}{3/8}\,dx=\frac83\left.\left(\frac{x^3}3- \frac{x^2}8\right)\right|_{-1}^1=\frac{16}9.$$ Мы же воспользуемся равенствами \eqref{3.3}. Имеем \begin{align*} \int_{-1}^1\,dx&=A_0+A_1+A_2,\\ \int_{-1}^1x\,dx&=A_0\left(-\frac12\right)+A_2\frac14,\\ \int_{-1}^1x^2\,dx&=A_0\left(-\frac12\right)^2+A_2\left(\frac14\right)^2. \end{align*} Таким образом, получаем следующую систему $$\begin{cases}3A_0+A_1+A_2&=2,\\ -\dfrac12A_0+\dfrac14A_2&=0,\\ \dfrac14A_0+\dfrac1{16}A_2&=\dfrac23.\end{cases}$$ Решая эту систему, находим $$A_0=\frac{16}9,\qquad A_1=-\frac{10}3,\qquad A_2=\frac{32}9.$$ Следовательно, искомая квадратурная формула имеет вид $$\int_{-1}^1f(x)\,dx\approx \frac{16}9f(-1/2)-\frac{10}3f(0)+\frac{32}9f(1/4 ).$$ \section{Квадратурные формулы Ньютона--Котеса} Для приближенного вычисления интеграла $$\int_a^bf(x)\,dx$$ разобьем отрезок $[a,b]$ на $n$ равных частей точками $$x_i=x_0+ih,\quad i=0,1,\ldots,n,\quad h=\frac{b-a}n$$ ($x_0=a$, $x_n=b$). Заменим функцию $f(x)$ на ее интерполяционный многочлен Лагранжа, построенный по значениям этой функции в точках $x_0,x_1,\ldots,x_n$. Для равноотстоящих узлов интерполяции (см.\ [7, стр.\ 12]) $$l_{ni}(x)=\frac{(-1)^{n-i}}{i!(n-i)!}q(q-1)\dots(q-i+1)(q-i-1)\dots(q-n),$$ где \begin{equation}\label{4.1} q=\frac{x-x_0}h. \end{equation} Сделав замену переменной \eqref{4.1} и учитывая, что $dx=h\,dq$, получим квадратурную формулу \eqref{3.1}, в которой $$A_i=h\frac{(-1)^{n-i}}{i!(n-i)!}\int_0^nq(q-1)\ldots(q-i+1)(q-i-1)\dots(q-n) \,dq.$$ Положим \begin{equation}\label{4.2} H_i=\frac1n\frac{(-1)^{n-i}}{i!(n-i)!}\int_0^nq(q-1)\ldots(q-i+1)(q-i-1) \dots(q-n)\,dq. \end{equation} Так как $h=\dfrac{b-a}n$, то $A_i=(b-a)H_i$. Тем самым получаем квадратурную формулу \begin{equation}\label{4.3} \int_a^bf(x)\,dx\approx(b-a)\sum_{i=0}^nH_if_i, \end{equation} называемую {\it квадратурной формулой Ньютона--Котеса}. Величины \eqref{4.2} называются {\it коэффициентами Котеса}. Для них существуют таблицы. Отметим некоторые простейшие свойства коэффициентов Котеса. Поскольку квадратурная формула \eqref{4.3} точна для любого многочлена степени $k\le n$, то она, в частности точна для $f(x)\equiv1$. Подставив эту функцию в \eqref{4.3}, находим \begin{equation}\label{4.4} \sum_{i=0}^nH_i=1. \end{equation} Кроме того, непосредственно из \eqref{4.2} вытекают равенства \begin{equation}\label{4.5} H_{n-i}=H_i,\quad i=0,1,\ldots,n-1. \end{equation} \section{Формула трапеций} Рассмотрим квадратурную формулу Ньютона--Котеса при $n=1$. В этом случае $h=b-a$, а сама формула имеет вид $$\int_a^bf(x)\,dx\approx h(H_0f_0+H_1f_1).$$ Пользуясь свойствами \eqref{4.4} и \eqref{4.5}, находим \begin{gather*} H_0+H_1=1,\\ H_0=H_1. \end{gather*} Следовательно, $H_0=H_1=1/2$. Итак, получена формула $$\int_a^bf(x)\,dx\approx h\frac{f_0+f_1}2,$$ называемая {\it формулой трапеций}. Оценим погрешность формулы трапеций. Для погрешности приближения функции многочленом Лагранжа первой степени имеет место равенство (см.\ [7, стр.\ 10, 12]) $$f(x)-L_1(x)=\frac{f''(\xi)}2h^2q(q-1),$$ где $\xi$ некоторая точка из интервала $(a,b)$, а $q$ определено равенством \eqref{4.1}. Таким образом, имеем \begin{multline*} \left|\int_a^bf(x)\,dx-h\frac{f_0+f_1}2\right|=\left|\int_a^b (f(x)-L_1(x))\,dx\right|\\ =\left|\int_0^1\frac{f''(\xi)}2h^3q(q-1)\,dq\right|\le\frac{M_2}2h^3\int_0^1q (1-q)\,dq=\frac{M_2}{12}h^3. \end{multline*} Аналогично усложненной формуле прямоугольников можно построить усложненную формулу трапеций, разбив весь отрезок $[a,b]$ на $n$ равных частей и применив на каждой из частей формулу трапеций. Тогда $h=\dfrac{b-a}n$, а соответствующая квадратурная формула имеет вид $$\int_a^bf(x)\,dx\approx h\left(\frac{f_0}2+f_1+\ldots+f_{n-1}+\frac{f_n}2 \right).$$ При этом для погрешности этой квадратурной формулы (мы обозначаем ее через $R_t$) будут справедливы неравенства \begin{equation}\label{5.1} R_t\le n\frac{M_2}{12}h^3=\frac{M_2}{12}(b-a)h^2=\frac{M_2(b-a)^3}{12n^2}. \end{equation} \begin{example}\label{p5.1}\rm В скольких точках надо вычислить функцию $f(x)=e^{x^2}$, чтобы при вычислении интеграла $$\int_0^1e^{x^2}\,dx$$ по усложненной формуле трапеций погрешность не превосходила $10^{-5}$. \end{example} {\sc Решение.} Поскольку $f''(x)=(2+4x^2)e^{x^2}$, то $M_2=6e$. Из \eqref{5.1} вытекает, что для достижения требуемой точности достаточно выполнения неравенства $$\frac{M_2}{12n^2}\le10^{-5}.$$ Отсюда $n\ge100\sqrt{5e}=368,6\ldots\,\,$. Число вычислений функции при разбиении отрезка на $n$ частей равно $n+1$, поэтому для вычисления рассматриваемого интеграла с заданной точностью потребуется находить значения функции в 370 точках. \section{Формула Симпсона} Пусть теперь $n=2$, т.е.\ $x_0=a$, $x_1=a+h$, $x_2=b$, а $h=\dfrac{b-a}2$. Из \eqref{4.2} имеем $$H_0=\frac14\int_0^2(q-1)(q-2)\,dq=\frac14\left.\left(\frac{q^3}3-\frac{3q^2 }2+2q\right)\right|_0^2=\frac16.$$ Тогда в силу \eqref{4.5} $H_2=H_0=1/6$, а из \eqref{4.4} вытекает, что $H_1=2/3$. Итак, мы получили квадратурную формулу $$\int_a^bf(x)\,dx\approx\frac h3(f_0+4f_1+f_2),$$ которая называется {\it формулой Симпсона}. Геометрический смысл этой формулы заключается в том, что вместо интеграла от исходной функции вычисляется интеграл от параболы, проходящей через значения функции в точках $a$, $\dfrac{a+b}2$ и $b$. Оценим погрешность формулы Симпсона. Без ограничения общности можно считать, что начало координат выбрано в точке $\dfrac{a+b}2$. Тем самым будем оценивать квадратурную формулу \begin{equation}\label{6.1} \int_{-h}^hf(x)\,dx\approx\frac h3(f(-h)+4f(0)+f(h)). \end{equation} Прежде всего отметим, что формула Симпсона точна не только на многочленах второй степени, как следует из замечания в п.\ 3, но и на многочленах третьей степени. Для того, чтобы убедиться в этом, достаточно проверить, что формула \eqref{6.1} точна для $f(x)=x^3$. Пусть $c$ --- произвольная точка из интервала $(0,h)$. Рассмотрим интерполяционный многочлен Лагранжа $L_3(x)$, интерполирующий функцию $f(x)$ в точках $-h$, $0$, $c$ и $h$. Имеем $$\int_{-h}^hL_3(x)\,dx=\frac h3(L_3(-h)+4L_3(0)+L_3(h))=\frac h3(f(-h)+4f(0) +f(h)).$$ Поэтому для погрешности формулы Симпсона, которую мы обозначим через $R$, будем иметь \begin{multline*} R=\biggl|\int_{-h}^hf(x)\,dx-\frac h3(f(-h)+4f(0)+f(h)) \biggr|\\ =\biggl|\int_{-h}^hf(x)\,dx-\int_{-h}^hL_3(x)\,dx\biggr|=\biggl|\int_{-h}^h(f (x)-L_3(x))\,dx\biggr|\\ \le\int_{-h}^h|f(x)-L_3(x)|\,dx. \end{multline*} Положим $$M_4=\max_{x\in[-h,h]}|f^{(4)}(x)|.$$ Тогда для погрешности приближения функции $f(x)$ многочленом Лагранжа справедлива следующая оценка (см.\ [7, стр.\ 11]) $$|f(x)-L_3(x)|\le\frac{M_4}{4!}|(x+h)x(x-c)(x-h)|.$$ Тем самым для любого $c\in(0,h)$ $$R\le\frac{M_4}{24}\int_{-h}^h|(x+h)x(x-c)(x-h)|\,dx.$$ Переходя к пределу в этом неравенстве при $c\to0$, получаем \begin{equation}\label{6.2} R\le\frac{M_4}{24}\int_{-h}^hx^2(h^2-x^2)\,dx=\frac{M_4}{24}\left.\left(h^2 \frac{x^3}3-\frac{x^5}5\right)\right|_{-h}^h=\frac{M_4}{90}h^5. \end{equation} Выведем теперь усложненную формулу Симпсона. Разделим отрезок $[a,b]$ на четное число отрезков $2n$ $$x_i=a+ih,\quad i=0,1,\ldots,2n,\quad h=\frac{b-a}{2n}.$$ На каждом из отрезков $[x_{2i},x_{2i+2}]$ применим формулу Симпсона $$\int_{x_{2i}}^{x_{2i+2}}f(x)\,dx\approx\frac h3(f_{2i}+4f_{2i+1}+f_{2i+2}).$$ Сложив эти равенства, получим усложненную формулу Симпсона $$\int_a^bf(x)\,dx\approx\frac h3\biggl(f_0+f_{2n}+4\sum_{i=1}^nf_{2i-1}+2 \sum_{i=1}^{n-1}f_{2i}\biggr).$$ Учитывая оценку \eqref{6.2}, для погрешности усложненной формулы Симпсона, которую мы обозначим через $R_s$, будем иметь $$R_s\le n\frac{M_4}{90}h^5=\frac{M_4}{180}(b-a)h^4,$$ где $$M_4=\max_{x\in[a,b]}|f^{(4)}(x)|.$$ Поскольку $h=\dfrac{b-a}{2n}$, то оценка величины $R_s$ может быть дана через $n$ $$R_s\le\frac{M_4(b-a)^5}{2880n^4}.$$ \begin{example}\rm Рассмотреть пример~\ref{p5.1} для формулы Симпсона. \end{example} {\sc Решение.} Имеем $f^{(4)}(x)=(12+48x^2+16x^4)e^{x^2}$. Тем самым $M_4=76e$. Для достижения требуемой точности достаточно выполнения неравенства $$\frac{M_4}{2880n^4}\le10^{-5}.$$ Отсюда $$n\ge10\sqrt[4]{\frac{19e}{72}}=9,2\ldots\,\,.$$ При использовании формулы Симпсона отрезок разбивается на $2n$ частей. Поэтому для вычисления рассматриваемого интеграла с заданной точностью потребуется находить значения функции в 21 точке (напомним, что в аналогичном примере для формулы трапеций при той же точности требовались вычисления в 370 точках). \section{Главная часть погрешности квадратурных формул} Для усложненной формулы прямоугольников было доказано следующее равенство \begin{equation}\label{7.1} \int_a^bf(x)\,dx=h\sum_{i=0}^{n-1}f_{i+1/2}+r_1, \end{equation} где $h=\dfrac{b-a}n$, а $$|r_1|\le\frac{M_2}{24}(b-a)h^2.$$ Мы получим более точную оценку в предположении, что функция $f(x)$ имеет непрерывную четвертую производную. Рассмотрим сначала случай, когда $[a,b]=[-h/2,h/2]$. По формуле Тейлора \begin{equation}\label{7.2} f(x)=f_0+f'_0x+\frac{f_0''}{2!}x^2+\frac{f_0'''}{3!}x^3+\frac{f^{(4)}(\xi)} {4!}x^4, \end{equation} где $f^{(i)}_0=f^{(i)}(0)$, $i=0,1,2,3$, а $\xi$ --- некоторая точка из интервала $(-h/2,h/2)$. Интегрируя \eqref{7.2}, получаем равенство $$\int_{-h/2}^{h/2}f(x)\,dx=f_0h+\frac{f_0''}{24}h^3+r_2,$$ в котором $$r_2=\int_{-h/2}^{h/2}\frac{f^{(4)}(\xi)}{4!}x^4\,dx.$$ Для величины $r_2$ справедлива следующая оценка $$|r_2|\le\frac{M_4}{4!}\int_{-h/2}^{h/2}x^4\,dx=\frac{M_4}{1920}h^5.$$ Таким образом, для усложненной формулы прямоугольников имеем \begin{equation}\label{7.3} \int_a^bf(x)\,dx=h\sum_{i=0}^{n-1}f_{i+1/2}+\frac{h^3}{24}\sum_{i=0}^{n-1}f ''_{i+1/2}+r_3, \end{equation} где $$|r_3|\le n\frac{M_4}{1920}h^5=\frac{M_4}{1920}(b-a)h^4.$$ Применяя \eqref{7.1} для второй производной, находим $$\int_a^bf''(x)\,dx=h\sum_{i=0}^{n-1}f''_{i+1/2}+r_4,$$ при этом $$|r_4|\le\frac{M_4}{24}(b-a)h^2.$$ Следовательно, $$h\sum_{i=0}^{n-1}f''_{i+1/2}=\int_a^bf''(x)\,dx-r_4.$$ Подставляя это выражение в \eqref{7.3}, получаем $$\int_a^bf(x)\,dx=h\sum_{i=0}^{n-1}f_{i+1/2}+\frac{h^2}{24}\int_a^bf''(x)\,d x+r_5,$$ где $$|r_5|=\Bigl|-\frac{h^2}{24}r_4+r_3\Bigr|\le\frac{M_4}{24^2}(b-a)h^4+\frac{M _4}{1920}(b-a)h^4=\frac{13M_4}{5760}(b-a)h^4.$$ {\sc Определение.} Говорят, что $\varphi(h)=O(h^k)$ (читается ``о большое''), если существует такая постоянная $C>0$, для которой $$|\varphi(h)|\le Ch^k.$$ Тем самым мы доказали равенство \begin{equation}\label{7.4} \int_a^bf(x)\,dx=h\sum_{i=0}^{n-1}f_{i+1/2}+ch^2+O(h^4), \end{equation} в котором $$c=\frac1{24}\int_a^bf''(x)\,dx.$$ Положим $$I=\int_a^bf(x)\,dx,\quad I_h^r=h\sum_{i=0}^{n-1}f_{i+1/2}.$$ Тогда равенство \eqref{7.4} запишется в виде $$I=I_h^r+ch^2+O(h^4).$$ Величина $ch^2$ называется {\it главной частью\/} погрешности формулы прямоугольников. Аналогичные равенства можно получить для формул трапеций и Симпсона \begin{align*} I&=I_h^t+c_1h^2+O(h^4),\\ I&=I_h^s+c_2h^4+O(h^6) \end{align*} (в последнем случае надо требовать, чтобы функция $f(x)$ имела непрерывную шестую производную). \section{Правило Рунге практической оценки погрешности} Пусть $z$ --- неизвестное точное значение некоторой величины, $z_h$ --- известное ее приближенное значение, зависящее от положительного параметра $h$, который может принимать сколь угодно малые значения. Предположим, что установлена связь между точным и приближенным значениями \begin{equation}\label{8.1} z=z_h+ch^k+O(h^{k+m}), \end{equation} где $c$ --- неизвестная не зависящая от $h$ постоянная. Тогда \begin{equation}\label{8.2} z=z_{h/2}+c\left(\frac h2\right)^k+O(h^{k+m}), \end{equation} так как для любой постоянной $C$ \ $O((Ch)^n)=O(h^n)$. Вычитая из \eqref{8.1} равенство \eqref{8.2}, будем иметь \begin{equation}\label{8.3} z_{h/2}-z_h=c\left(\frac h2\right)^k(2^k-1)+O(h^{k+m}). \end{equation} Отсюда $$c\left(\frac h2\right)^k=\frac{z_{h/2}-z_h}{2^k-1}+O(h^{k+m}).$$ Следовательно, при $c\ne0$ величина $\dfrac{z_{h/2}-z_h}{2^k-1}$ отличается от главного члена погрешности $z-z_{h/2}$ на $O(h^{k+m})$, т.е. \begin{equation}\label{8.4} z-z_{h/2}=\frac{z_{h/2}-z_h}{2^k-1}+O(h^{k+m}). \end{equation} Тем самым при $c\ne0$ оценить погрешность можно так: $$z-z_{h/2}\approx\frac{z_{h/2}-z_h}{2^k-1}.$$ Такой способ оценки погрешности называется {\it правилом Рунге}. {\sc Замечание.} На практике считается, что условие $c\ne0$ выполнено, если \begin{equation}\label{8.5} \left|2^k\frac{z_{h/2}-z_h}{z_h-z_{2h}}-1\right|<0,1. \end{equation} Только в этом случае рекомендуется применение правила Рунге. Поясним условие \eqref{8.5}. Из \eqref{8.3} следует, что \begin{align*} z_{h/2}-z_h&=c\left(\frac h2\right)^k(2^k-1)+O(h^{k+m}),\\ z_h-z_{2h}&=ch^k(2^k-1)+O(h^{k+m}). \end{align*} Поэтому $$\frac{z_{h/2}-z_h}{z_h-z_{2h}}\approx\frac1{2^k}.$$ Таким образом, $$2^k\frac{z_{h/2}-z_h}{z_h-z_{2h}}\approx1.$$ \section{Уточнение приближенного решения по Ричардсону} Положим $$z_h^*=\frac{2^kz_{h/2}-z_h}{2^k-1}.$$ Тогда из \eqref{8.4} получаем $$z=z_{h/2}+\frac{z_{h/2}-z_h}{2^k-1}+O(h^{k+m})=z_h^*+O(h^{k+m}).$$ При $c\ne0$ \ $z-z_{h/2}$ имеет $k$-ый порядок малости, а $z-z_h^*$ --- $(k+m)$-ый порядок малости, т.е.\ $z_h^*$ --- более точное приближение. Оно носит название {\it уточнение по Ричардсону}. Таким образом, если вычисляется определенный интеграл по усложненным формулам прямоугольников или трапеций, то $k=2$, и мы можем оценить приближенно погрешность по правилу Рунге $$I-I_{h/2}\approx\frac{I_{h/2}-I_h}3.$$ Кроме того, можно найти уточнение по Ричардсону $$I_h^*=\frac{4I_{h/2}-I_h}3.$$ Для формулы Симпсона $k=4$, и приближенная оценка погрешности имеет вид $$I-I^s_{h/2}\approx\frac{I_{h/2}^s-I_h^s}{15}.$$ Для уточнения по Ричардсону имеем $${I_h^s}^*=\frac{16I_{h/2}^s-I_h^s}{15}.$$ {\sc Замечание.} Формулы трапеций и Симпсона удобны тем, что при переходе от $h$ к $h/2$ все вычисленные ранее значения функций используются в новой квадратурной формуле. Правило Рунге и уточнение по Ричардсону можно применять и для других задач приближенного вычисления. Рассмотрим в качестве примера численное дифференцирование. Пусть известны значения некоторой достаточно гладкой функции $f(x)$ в точках $x_i=x_0+ih$, $i=0,\pm1$. Положим $f_i=f(x_i)$, $i=0,\pm1$. Хорошо известна формула численного дифференцирования (см., например, [7, стр.\ 26]) $$f_0''\approx\frac{f_1-2f_0+f_{-1}}{h^2}=f''_{0h}.$$ Найдем главную часть погрешности в этом методе. По формуле Тейлора имеем $$f_{\pm1}=f_0\pm f_0'h+\frac{f_0''}{2!}h^2\pm\frac{f_0'''}{3!}h^3+\frac{f_0^ {(4)}}{4!}h^4\pm\frac{f_0^{(5)}}{5!}+\frac{f^{(6)}(\xi_{\pm1})}{6!}h^6,$$ где $f_0^{(i)}=f^{(i)}(0)$, $i=0,1,\ldots,5$, $\xi_1\in(0,h)$, а $\xi_{-1}\in(-h,0)$. Положим $$r=-\frac{f^{(6)}(\xi_1)+f^{(6)}(\xi_{-1})}{6!}h^4.$$ Тогда $$f_1+f_{-1}=2f_0+f_0''h^2+\frac{f_0^{(4)}}{12}h^4-rh^2.$$ Тем самым \begin{equation}\label{9.1} f_0''=\frac{f_1-2f_0+f_{-1}}{h^2}-\frac{f_0^{(4)}}{12}h^2+r. \end{equation} Поскольку $$|r|\le\frac{2M_6}{6!}h^4,$$ то равенство \eqref{9.1} может быть записано в виде $$f_0''=f''_{0h}+ch^2+O(h^4),$$ где $c=-f_0^{(4)}/12$. Следовательно, по правилу Рунге $$f_0''-f''_{0,h/2}\approx\frac{f''_{0,h/2}-f''_{0h}}3,$$ а для уточнения по Ричардсону справедливо равенство $${f''_{0h}}^*=\frac{4f''_{0,h/2}-f''_{0h}}3.$$ \section*{\bf Литература} \begin{itemize} \item[1.] Бахвалов Н.С. {\it Численные методы}. М.: Наука, 1973. \item[2.] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. {\it Численные методы}. М.: Наука, 1987. \item[3.] Березин И.С., Жидков Н.П. {\it Методы вычислений}. М.: Наука, 1966. Т.1; Физматгиз, 1962. Т.2. \item[4.] Волков Е.А. {\it Численные методы}. М.: Наука, 1982. \item[5.] Демидович Б.П., Марон И.А. {\it Основы вычислительной математики}. М.: Наука, 1966. \item[6.] Калиткин Н.Н. {\it Численные методы}. М.: Наука, 1978. \item[7.] Осипенко К.Ю. {\it Аппроксимация функций многочленами и численное дифференцирование}: Методические указания по курсу ``Численные методы''; МГАТУ. М., 1994. \end{itemize} \tableofcontents \end{document}