%\documentstyle{report}
%\documentclass{report }

%\usepackage[cp866]{inputenc}
%\usepackage[cp1251]{inputenc}

%\usepackage[russian]{babel}
%\usepackage[LCY]{fontenc}
%\usepackage{graphicx,mybib,russian}
%\newcommand{\<<}{\lfq}
\documentclass[12pt,a4paper]{article}
%\usepackage[cp866]{inputenc}
\usepackage{pscyr}
\usepackage[cp1251]{inputenc}

\usepackage[russian]{babel}
%\usepackage[LCY]{fontenc}
%\usepackage{graphicx,mybib,russian}
%\newcommand{\<<}{\lfq}
\textheight 22cm \textwidth 16cm \topmargin 0cm \oddsidemargin 0cm \evensidemargin 0cm
%\twocolumn
%\oddsidemargin =-2cm
%\topmargin -3cm
%\textheight=27cm
%\textwidth=17cm
\pagestyle{empty}
\begin{document}

{\footnotesize

Лекция 11.{\bf Усложненные квадратурные формулы и их погрешность}

Мы видели, что в квадратурных формулах Ньютона-Котеса нельзя для уменьшения
их погрешности  увеличивать степень интерполяционного многочлена.
Оказывается, что в этом и нет необходимости т.к. погрешность можно уменьшить за счет
другого приема - предварительного разбиения интервала интегрирования $[a,b]$ и
построенные таким образом квадратурные формулы носят название усложненных квадратурных
формул.

Для построения усложненных квадратурных формул
для интегрирования функций $f(x)\in C{[a,b]}^{n+1}$ с весовой
$p(x)$-заданного вида, построим следующую процедуру

$ \int_a^bf(x)p(x)dx=\{[a,b]\to \omega_n=\{x_i^{(n)}:a\le x_0\le
x_1\cdots\le x_n\le b \} \mbox{и интеграл представим в виде} \}=
\sum_{i=0}^{n-1}\int_{x_i}^{x_{i+1}}f(x)p(x)dx=
\{ \mbox{ замена переменных}
$

$
 x=x_i+(x_{i+1}-x_i)t \}=
\sum_{i=0}^{n-1}(x_{i+1}-x_i)\int_0^1f(x_i+(x_{i+1}-x_i)t)p(x_i+(x_{i+1}-x_i)t)dt =
\{[0,1]\to \omega_M=\{t_i:t_0\le 0\le \cdots\le t_M\le 1   \}$

$
\mbox{на этой сетке построим квадратурную}
$

$
\mbox{ формулу Ньютона-Котеса } \}\to
$

$
\sum_{i=0}^{n-1}(x_{i+1}-x_i)\sum_{k=0}^{M}f(x_i+(x_{i+1}-x_i)t_k)$

$
\int_0^1\frac{\omega_k (t)}{t-t_k}\omega '(t)p(x_i+(x_{i+1}-x_i)t)dt+
\sum_{i=0}^{n-1}(x_{i+1}-x_i)\sum_{k=0}^M\int_0^1
\frac{f^{M+1}(\xi)}{(M+1)!}p(x_i+(x_{i+1}-x_i)t)
\omega_k(t)dt, \xi\in [0,1]
$

Первое слогаемое в этой формуле - квадратурная формула, второе -ее погрешность.


{\bf Простейшие усложненные квадратурные

формулы и их погрешность для $p(x)=1$}

1. Квадратурная формула трапеций.

$ M=1,\omega_M=\{0,1\}$

$\int_a^bf(x)dx\to\sum_{i=0}^{n-1}(x_{i+1}-x_i)\frac{1}{2}[f(x_i)+f(x_{i+1})]$

ее погрешность

$D_n^{tr}=\sum_{i=0}^{n-1}\frac{f'''(\xi)}{12}(x_{i+1}-x_i)^3$

$h=|x_{i+1}-x_i|$ -однородные сетки $hn=(b-a), \sum_{i=0}^{n-1}h=hn=b-a$

имеем

$\int_a^bf(x)dx\to h[\frac{1}{2}f(x_0)+f(x_1)+\ldots +f(x_{n-1})+\frac{1}{2}f(x_n)]$

$D_h^{tr}=\frac{b-a}{12}\max_{x\in [a,b]}|f'''(x)|h^2 ,f(x)\in C^2_{[a,b]}=O(h^2)$

Более точные оценки для функций $f(x)\in C^4_{[a,b]}$ дают для порядка малости

$D_h^{tr}=O(h^2)+O(h^4)$.

2. Усложненная квадратурная формула Симпсона.

$M=2,\omega_M=\{0,1/2,1\}$

$\int_a^bf(x)dx\to\sum_{i=0}^{n-1}(x_{i+1}-x_i)\frac{1}{6}[f(x_i)+
4*f(\frac{x_i+x_{i+1}}{2})+f(x_{i+1})]$

ее погрешность

$D_n^{Sim}=\sum_{i=0}^{n-1}\frac{f^{(5)}(\xi)}{2880}(x_{i+1}-x_i)^5$

$h=|x_{i+1}-x_i|$ -однородные сетки $hn=(b-a), \sum_{i=0}^{n-1}h=hn=b-a$
$n\to 2n$ -четные
имеем
$\int_a^b f(x)dx\to \frac{h}{3}[f(x_0)+4f(x_1)+
2f(x_2) \ldots +4f(x_{2n-1})+f(x_{2n})] $

$D_h^{Sim}=\frac{b-a}{2880}\max_{x\in [a,b]}|f^{(4)}(x)|h^4 ,f(x)\in C^4_{[a,b]}=O(h^2)$

Более точные оценки для функций $f(x)\in C^6_{[a,b]}$ дают для порядка малости

$D_h^{Sim}=O(h^4)+O(h^6)$.


3. Сплайн- квадратура.
 Пусть $f(x)\in C^4_{[a,b]}$

$\int_a^bf(x)dx=\{ [a,b]\to \omega_n=\{x_i: a\le x_0<x_1\ldots<x_n\le b   \}$
{сеточную функцию заменим сплайном 3-степени
}$S_3=a_i+b_i(x-x_{i-1})+c_i(x-x_{i-1})^2+d_i(x-x_{i-1})^3 \}\to
\sum_{i=0}^{n-1}[a_ih_i+b_i\frac{h_i^2}{2}+c_i\frac{h_i^3}{3}+d_i\frac{h_i^4}{4}]+
O(h^4)$

$h_i=|x_i-x_{i-1}|, h=\max_{i}h_i$.

Если вспомнить определение кубического сплайна, то

$a_i=f(x_i), b_i=\frac{f(x_i)-f(x_{i-1})}{h_i}-\frac{h_i}{3}(c_{i+1}+2c_i),
d_i=\frac{c_{i+1}-c_i}{3h_i}$,  коэффициенты $c_i$ -находятся из системы
линейных уравнений с матрицей ленточного типа, методом прогонки

$\int_a^bf(x)dx\to \sum_{i=0}^{n-1}[\frac{h_i}{2}(f(x_i)+f(x_{i+1})-
\frac{h_i^3}{3}(c_i+c_{i+1})+O(h^4)$

{\bf Правило Рунге и уточнение решения по Ричардсону}.

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

Пусть на сетке $\omega_n=\{x_i:,a\le x_0<x_1\cdots<x_n\le b\}, h=\max_i|x_{i+1}-x_i |<1$
известно, что точное значение $I$ некоторой величины связано с ее приближенным
значением $I_h$ связью вида

$I=I_h+\alpha_qh^q+\alpha_ph^p , p>q, \alpha_q,\alpha_p$-посттоянные, не зависящие от $h$.

$(1) \omega_{h_1} \to I=I_{h_1}+\alpha_q h_1^q+\alpha_p h_1^p=I_{h_1}+D_{h_1} $


$(2) \omega_{h_2}\to I=I_{h_2}+\alpha_qh_2^q+\alpha_ph_2^p=I_{h_2}+D_{h_2} $

Вычитая из (1) (2), будем иметь

$I_{h_1}-I_{h_2}=\alpha_qh_2^q[(\frac{h_1}{h_2})^q-1]+
\alpha_ph_2^p[(\frac{h_1}{h_2})^p-1]$

Откуда следует, что

$\alpha_qh^q_2=\frac{I_{h_1}-I_{h_2}-\alpha_ph_2^p[(\frac{h_1}{h_2})^p-1]
}{[(\frac{h_1}{h_2})^q-1]}\to (2)$  будем иметь

$I-I_{h_2}=\frac{I_{h_1}-I_{h_2}-\alpha_ph_2^p[(\frac{h_1}{h_2})^p-1]
}{[(\frac{h_1}{h_2})^q-1]}+ \alpha_ph_2^p\sim
\frac{I_{h_1}-I_{h_2}}{[(\frac{h_1}{h_2})^q-1]}+O(h^P)$
Откуда следует {\bf правило}, если на сетке $\omega_h$ точное значение некоторой величины
$I$ связано с приближенным значением $I_h$ соотношением
$I=I_h+\alpha_qh^q+\alpha_ph^p$ и приближенные значения $I_{h_1},I_{h_2}$, вычисленные на сетках
$\omega_{h_1},\omega_{h_2}$ удовлетворяют соотношению

$\frac {| I_{h_1}-I_{h_2} | }{ |(\frac{h_1}{h_2})^q-1|}
<\epsilon $, $\epsilon $-наперед заданное
малое число, тогда

$ |I-I_{h_2}|<\epsilon $.

Кроме того, оказывается, что с помощью линейной комбинации этих же приближенных значений,
 можно получить приближенное значение с порядком малости $O(h^p)$. Действительно

$I=\gamma I+(1-\gamma)I=\{(1),(2)\to\}=\gamma I_{h_1}+(1-\gamma)I_{h_2}+
\gamma D_{h_1}+(1-\gamma)D_{h_2}$

$\gamma D_{h_1}+(1-\gamma)D_{h_2}=\gamma \alpha_qh_1^q+(1-\gamma)\alpha_qh_2^q+
\gamma \alpha_ph_1^p+(1-\gamma)\alpha_ph_2^p$

Приравнивая нулю слагаемое порядка малости $O(h^q)$

$\gamma \alpha_qh_1^q+(1-\gamma)\alpha_qh_2^q=0$ находим, что

$\bar{\gamma}=\frac{h_2^q}{h_2^q-h_1^q}$  и, следовательно величина

$I_h=\bar{\gamma}I_{h_1}+(1-\bar{\gamma})I_{h_2}$ будет отличатся от точного
значения на величину $O(h^p)$. Найденное таким образом приближенное значение и
называется уточнением по Ричардсону.

На практике $h_2=h_1/2,\qquad \gamma=0.5.$

 Правило Рунге на практике следует применять с осторожностью, т.к. например
в процессе вычисления интегралов могу встречаться функции с дельтообразной
особенностью для которых выбор сеток неочевиден. (графический рис.)
В стандартных программах, когда говорят о вычислении интеграла с заданной точностью
имеют в виду применение правила Рунге и соответствующие алгоритмы называют
адаптивными.

{\bf Интегралы с особенностями} Очевидно, что несобственные интегралы, не вычисляются
по простейшим квадратурным формулам и для преодоления трудностей их вычисления часто помогают
различные приемы, например, аналитические методы

1. Интегрирование по частям

$\int_0^{10}\frac{\cos(x)}{\sqrt{x}}dx=\cos(x)(2\sqrt{x})|_0^{10}+
2\int_0^{10}\sqrt(x)\sin(x)dx$

в данном случае избавились от особенности

2. Замена переменной

$\int_0^c\frac{\cos(x)}{\sqrt(x)}dx=\{x=t^2\}=\int_0^{c^2}\frac{\cos(t^2)}{t}2tdt=
$

$
2\int_0^{c^2}\cos(t^2)dt$

 в данном случае мы избавились от особенности,
но получили неприятности при вычислении интеграла для больших $c$.

3. Использование рядов.

$\int_0^{10}\frac{\cos(x)}{\sqrt{x}}dx=
\int_0^{c}\frac{\cos(x)}{\sqrt{x}}dx+\int_c^{10}\frac{\cos(x)}{\sqrt{x}}dx$

И для вычисления первого интеграла с особенностью

$\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}+\cdots$

Тогда

$ \int_0^c \frac{\cos(x)}{\sqrt{x}}dx=\int_0^c\frac{1}{\sqrt{x}}dx-
\int_0^c \frac{x^{3/2}} {2!}dx +\int_0^c \frac{x^{7/2}}{4!}dx-\cdots $

Каждый из интегралов в правой части вычисляется по формуле Ньютона-Лейбница и
вопрос вычисления интеграла с особенностью сводится к вычислению ряда.


{\bf Использование весовых функций}

Очень часто для вычисления интегралов требуется выделить весовую функцию
и построить квадратурную формулу в соответствующем виде. Например, интегралы
от быстро осциллирующих функций

$\int_a^bf(x)\sin(\omega x)dx, \int_a^bf(x)\cos(\omega x)dx$

Здесь удобно интерполяционным многочленом заменить функцию $f(x)$, а
функции $\sin(\omega x),\cos(\omega x)$ взять в качестве весовых.


{\bf Кратные интегралы}
Вычисление интеграла по некоторой области $D$

$I^{(m)}=I^{(m)}(f,D)=\int\cdots\int\limits_Df(x_1,x_2,
\cdots,x_m)dx_1\cdots dx_m$

можно свести  к вычислению определенных интегралов $\int_a^bf(x)dx$
Однако сложность соответствующих алгоритмов и число арифметических
операций для вычисления интеграла с заданной точностью возрастает с ростом размерности
и дополнительное затруднение возникает для областей общего вида.
Наиболее распространенным в этом случае алгоритме является метод Монте-Карло.
Суть метода в следующем. Будем считать, что область интегрирования
лежит внутри единичного куба $D_1=\{|x_i|< 1,i=\overline{1,n}\}$. Доопределим
функцию $f(x),x=(x_1<\cdots,x_n)$ вне области $D$ нулем. Пусть далее имеется
датчик случайных чисел, который выдает случайные числа, лежащие на отрезке
$[0,1]$ с плотностью распределения вероятности, равной единице. С помощью этого
датчика построим последовательность точек $P_k, k=1,2,\cdots $ и составим
числовую последовательность

$ S_k=\frac{1}{k}\sum_{j=0}^kf(P_j).$

Легко представить, что при произвольном $\epsilon>0 $ вероятность того, что
погрешность $I^{(m)}-S_k$ меньше $\epsilon$

$|I^{(m)}-S_k|<\epsilon$

$P(|I^{(m)}-S_k|,\epsilon)\to_{k\to 0} 1$

 с ростом числа $k$ испытаний стремится
к единице. Имеются соответствующие теоремы.

Погрешность вычислений при этом имеет порядок малости $O(\frac{1}{\sqrt{k}}).$


Смысл метода особенно нагляден если $n=2,f(x)=1$. Интеграл в этом случае -площадь
области $D$, а $S_k$ - отношение числа тех точек $P_1,P_2,\cdots,P_k$, котрые
попали в область $D$, к числу всех точек. Очевидно, что это отношение приблизительно равно
площади области $D$, поскольку площадь всего квадрата куда попадают точки, равна единице.

Привлекательность метода -в логической простоте алгоритма, который не усложняется
с ростом размерности и усложнением области интегрирования.
Однако в случае простых областей и узких классов функций метод проигрывает методам
учитывающих эту специфику.

}
\end{document}
