Численное интегрирование на разностных схемах

Опишем методику численного расчета интеграла \(\int_a^b f(x) d x\) с помощью разностных схем в старших порядках. В отличие от квадратур, тут нам будут нужны производные.

Теорема. Пусть функция \(f(x)\) аналитична в каждой точке отрезка \([a ; b]\). Тогда её интеграл можно вычислять квадратурой на разностных схемах старшего порядка по формуле:

\[ \begin{gathered} \int_a^b f(x) d x=h \sum_{j=0}^{J-1} \sum_{k=-m}^m W_{m-k}^m y_{j+k}+O\left(h^{2 m+2}\right) \end{gathered} \]

где коэффициенты не зависят от функции и задаются рекурсивными формулами:

\[ W_k^m=\sum_{n=0}^m \frac{A_{k, 2 n}^{m-n}}{2^{2 n}(2 n+1) !} \]

\[ A_{k n}^m=\sum_{l=0}^m(-1)^{k-m} D_n^{(l)} C_{n+2 l}^{k-m+l} \]

\[ D_n^{(j+1)}=\sum_{k=0}^{n+2 j} \sum_{l=0}^j(-1)^k D_n^{(l)} C_{n+2 l}^{k-j+l} \frac{(n+2 j-2 k)^{n+2 j+2}}{2^{n+2 j+2}(n+2 j+2) !}, \quad D_n^{(0)} \equiv 1 \]

Примечание. Обратите внимание, что квадратура вылезает за сетку на отрезке \([a ; b]\)! Это специфика этого метода. Про бесконечные пределы смотреть примечание ниже.

Доказательство.

Пусть функция \(f(x)\) аналитична в каждой точке отрезка \([a ; b]\). Тогда она аналитична хотя бы немного за пределами отрезка \([a ; b]\), т. е. на некотором интервале \(\left(a^{\prime} ; b^{\prime}\right)\supset[a ; b]\). Это условие нам понадобится для того чтобы использовать только симметричные разностные схемы даже на границах отрезка.

Разобьем отрезок \([a ; b]\) на \(J\) одинаковых кусочков и на каждом кусочке разложим \(f(x)\) в ряд Тейлора в средней точке каждого кусочка, после чего вычислим интеграл на каждом кусочке. Имеем,

\[\begin{split} \begin{gathered} \int_a^b f(x) d x=\sum_{j=0}^{J-1} \int_{x_j-\Delta x}^{x_j+\Delta x} f(x) d x= \\ =\sum_{j=0}^{J-1} \sum_{n=0}^m \frac{f^{(2 n)}\left(x_j\right)}{(2 n) !} \int_{x_j-\Delta x}^{x_j+\Delta x}\left(x-x_j\right)^{2 n} d x+O\left(\Delta x^{2 m+2}\right)= \\ =\sum_{j=0}^{J-1} \sum_{n=0}^m \frac{f^{(2 n)}\left(x_j\right)}{(2 n+1) !} 2 \Delta x^{2 n+1}+O\left(\Delta x^{2 m+2}\right) . \end{gathered} \end{split}\]

Где \(\Delta x=(b-a) / 2 J=\cfrac{h}{2}\)- половина длины кусочка, \(J-\) число кусочков, \(x_j=a+(2 j+1) \Delta x-\) средняя точка \(j\) го кусочка, и учтено, что \(\sum_{j=0}^{J-1} O(\Delta x)=O(1)\) и что интеграл в симметричных пределах от нечетных степеней \(\left(x-x_j\right)^{2 n+1}\) равен нулю.

Итого получили метод численного расчёта интеграла произвольного порядка сходимости \(2m+2\), в котором нам нужны не только значения функции в узлах сетки, но и значения всех производных до \(2m\) включительно в этих же узлах с нужной точностью.

В следующем году в семинаре 4.1 будет выведена общая формула для симметричной разностной схемы произвольного порядка аппроксимации. Пока её вывод можно посмотреть в файле numer.pdf в доп. материалах. Сама формула:

\[ f^{(n)}\left(x_j\right)=\frac{1}{2^n h^n} \sum_{k=0}^{n+2 m} A_{k n}^m \cdot y_{j+n+2 m-2 k}+O\left(h^{2 m+2}\right) \]

\[ A_{k n}^m=\sum_{l=0}^m(-1)^{k-m} D_n^{(l)} C_{n+2 l}^{k-m+l} \]

\[ D_n^{(j+1)}=\sum_{k=0}^{n+2 j} \sum_{l=0}^j(-1)^k D_n^{(l)} C_{n+2 l}^{k-j+l} \frac{(n+2 j-2 k)^{n+2 j+2}}{2^{n+2 j+2}(n+2 j+2) !}, \quad D_n^{(0)} \equiv 1 \]

Отметим, что это разностная схема не сходится с ростом \(m\)! при фиксированном шаге. Поэтому рекомендуется брать \(m\) меньше, чем точность исходных данных. Для машинной точности (аналитические функции, 16 значящих знаков) хорошее значение \(m=7\).

Нам нужны только чётные производные и даже немаксимального порядка аппроксимации (т.к. для каждой из них производные домножаются на \(\Delta x = \frac{h}{2}\) в каких-то степенях). Внесём нужные модификации в формулу:

\[ f^{(2 n)}\left(x_j\right)=\frac{1}{(2 \Delta x)^{2 n}} \sum_{k=0}^{2 m} A_{k, 2 n}^{m-n} \cdot y_{j+m-k}+O\left(\Delta x^{2 m-2 n+2}\right) \]

Подставляя эту формулу в оценку интеграла, имеем:

\[ \int_a^b f(x) d x=2 \Delta x \sum_{j=0}^{J-1} \sum_{n=0}^m \sum_{k=0}^{2 m} \frac{A_{k, 2 n}^{m-n}}{2^{2 n}(2 n+1) !} \cdot y_{j+m-k}+O\left(\Delta x^{2 m+2}\right) \]

Вводя коэфициенты

\[ W_k^m=\sum_{n=0}^m \frac{A_{k, 2 n}^{m-n}}{2^{2 n}(2 n+1) !} \]

получаем квадратуру, основанную на разностных схемах на \(2m\) точках (заметье, что она вылезает за отрезок интегрирования!):

\[\begin{split} \begin{gathered} \int_a^b f(x) d x&=2 \Delta x \sum_{j=0}^{J-1} \sum_{k=-m}^m W_{m-k}^m y_{j+k}+O\left(\Delta x^{2 m+2}\right)= \\ &=h \sum_{j=0}^{J-1} \sum_{k=-m}^m W_{m-k}^m y_{j+k}+O\left(h^{2 m+2}\right) \end{gathered} \end{split}\]

Примечание. Важно заметить, что квадратура для интеграла «вылезает» за пределы области интегрирования - отрезка \([a ; b]\). Это - специфика применяемых здесь разностных схем в старших порядках. Но зато, интеграл вычисляется очень быстро -\(J+2 m\) обращений к функции дают ошибку \(O\left(\Delta x^{2 m+2}\right)\). Так, при \(m=7\) вычисление сетки значений \(\exp \left(x_j\right)\) для интерграла \(\int_a^b e^x d x\) машинная точность \(10^{-16}\) достигается при длине кусочка \(x_{j+1}-x_j=0.25\). Тогда как обычный метод Симпсона (с ошибкой \(O\left(\Delta x^4\right)\) ) для достижения машинной точности требует шаг \(\Delta x=0.0001\), что на три порядка меньше!

Более того, для данной квадратуры естественно можно применять алгоритм Рунге для оценки погрешности.

Примечание. Отметим, что если выход сетки \(\left\{x_j\right\}\) за пределы области интегрирования недопустим (неаналитичность или бесконечные пределы), можно, как один из вариантов, вместо интеграла \(\int_a^b f(x) d x\) вычислять интеграл \(\int_0^1 f(x(t)) \dot{x}(t) d t\) с заменами, например,

\[x(t)=a+(b-a)(1-\cos (\pi t)) / 2\]

или при \(b=+\infty\)

\[ x(t)=\frac{2 a}{(1+\cos (\pi t))} \]

Аналитичность функции \(f(x(t)) \dot{x}(t)\) в каждой точке отрезка \(t \in[0 ; 1]\) должна гарантироваться заменой в обязательном порядке. В противном случае нужно использовать другую замену.