Расчёт бесконечных сумм

Содержание

Расчёт бесконечных сумм#

Метод прямой сшивки#

С помощью только что выведенного квадратуры можно вычислять очень большие и даже бесконечные суммы.

Преобразуем итоговую формулу в прошлом пункте к следующему виду (с учётом $\sum_{k=0}^{2 m} W_k^m=1$):

\[ \begin{gathered} \int_a^b f(x) d x=2 \Delta x\left[\sum_{j=m}^{J-m-1} y_j+\sum_{j=-m}^{m-1} I_{j+m}^m y_j+\sum_{j=J-m}^{J+m-1}\left(1-I_{j+m-J}^m\right) y_j\right]+O\left(\Delta x^{2 m+2}\right) \end{gathered} \]

где опять введены новые коэффициенты

\[ I_l^m=\sum_{k=0}^l W_{2 m-k}^m \]

Выразим сумму через интеграл:

\[ \sum_{j=m}^{J+m-1} y_j=\frac{1}{2 \Delta x} \int_a^b f(x) d x-\sum_{j=-m}^{m-1} I_{j+m}^m y_j+\sum_{j=J-m}^{J+m-1} I_{j+m-J}^m y_j + O\left(\Delta x^{2 m+1}\right) \]

Напомним, что $y_j=f(a+(2 j+1) \Delta x), m \lesssim 7-\text { порядок схемы, а } b=a+2 J \Delta x$.

Мы видим, что каково бы ни было большим число $J$, сумма большого слагаемых ( $J$ штук) выражается через сумму малого числа слагаемых ( $4 m$ штук) плюс хвостовой интеграл. Иначе говоря, сделав лишь $4 m$ обращений к последовательности, можно перейти от суммы к интегралу.

Ошибку перехода $O\left(\Delta x^{2 m+1}\right)$ можно опустить, когда последовательность $\left\{y_j\right\}$ меняется с $j$ медленно, в смысле, функция $f(x)$ меняется медленно на масштабе $\Delta x$. Это следует из независимости ошибки от $J$.

Теперь изложим сам метод.

Пусть нужно вычислить сумму $\sum_{j=j_0}^{j_1-1} a_j$, где величина $j_1-j_0$ может быть велика или даже равна бесконечности, а последовательность $\left\{a_j\right\}$ изменяется медленно с $j$ при больших $j$. Более точно, при $j>j_m$ характерный масштаб изменения последовательности $\left\{a_j\right\}$ больше или равен некоторому $j_{\mathrm{s}} \gg 1$.

Введём функцию $f(2 j \Delta x)=a_j, a=\left(2 j_m+1\right) \Delta x$ и $J=j_1-j_m-1$, предполагая необходимые аналитические свойства функции $f(x)$ (грубо, предполагаем работу с аналитическим продолжением $f(x)$). Тогда связь обозначений будет иметь вид $y_j=$ $a_{j+j_m-1}$, при этом $b=\left(2 j_1-1\right) \Delta x$. После такой смены обозначений:

\[ \int_a^b f(x) \frac{d x}{2 \Delta x}=\int_{-1 / 2}^{J-1 / 2} y_j d j=\int_{j_m+1 / 2}^{j_1-1 / 2} a_j d j \]

где $a_j$ и $y_j$ мы рассматриваем, как функции непрерывного аргумента $j$. Правый интеграл таким образом по квадратуре:

\[\begin{split} \begin{gathered} \int_{j_m+1 / 2}^{j_1-1 / 2} a_j d j= \\ =\sum_{j=j_m+m+1}^{j_1-m-1} a_j+\sum_{j=j_m-m+1}^{j_m+m} I_{j+m-j_m-1}^m a_j+\sum_{j=j_1-m}^{j_1+m-1}\left(1-I_{j+m-j_1}^m\right) a_j+O\left(\Delta x^{2 m+2}\right) \end{gathered} \end{split}\]

Ошибка $O\left(\Delta x^{2 m+2}\right)$ мала как $O\left(\Delta x / x_{\mathrm{s}}\right)^{2 m+2}$, где $x_{\mathrm{s}}$ - характерный масштаб изменения функции $f(x)$ на отрезке $[a ; b]$, включая краевые точки. В обозначениях этого параграфа $x_{\mathrm{s}}=2 j_{\mathrm{s}} \Delta x$. Отсюда после преобразований получаем искомое выражение

\[\begin{split} \begin{gathered} \sum_{j=j_0}^{j_1-1} a_j=\int_{j_m+1 / 2}^{j_1-1 / 2} a_j d j+\sum_{j=j_0}^{j_m-m} a_j-\sum_{j=j_1}^{j_1+m-1} a_j+ \\ +\sum_{j=j_m-m+1}^{j_m+m}\left(1-I_{j+m-j_m-1}^m\right) a_j+\sum_{j=j_1-m}^{j_1+m-1} I_{j+m-j_1}^m a_j+O\left(1 / j_{\mathrm{s}}^{2 m+2}\right) . \end{gathered} \end{split}\]

для суммы большого числа слагаемых (порядка $j_1-j_0$ штук) через хвостовой интеграл и сумму малого числа слагаемых (порядка $m$ и $j_m-j_0$ штук).

Сумма ряда 1. В важном частном случае $j_0=1$ и $j_1=\infty$ формула приобретает вид

\[ \sum_{j=1}^{\infty} a_j=\int_{j_m+1 / 2}^{\infty} a_j d j+\sum_{j=1}^{j_m-m} a_j+\sum_{j=j_m-m+1}^{j_m+m}\left(1-I_{j+m-j_m-1}^m\right) a_j+O\left(1 / j_{\mathrm{s}}^{2 m+2}\right) \]

Сумма ряда 2. В другом важном частном случае $j_0=-\infty$ и $j_1=\infty$, выбирая $j_m=-\infty$, получаем обоснование перехода от суммы к интегралу

\[ \sum_{j=-\infty}^{\infty} a_j=\int_{-\infty}^{\infty} a_j d j+O\left(1 / j_{\mathrm{s}}^{2 m+2}\right) . \]

для медленно меняющейся последовательности $\left\{a_j\right\}_{j=-\infty}^{\infty}$, характерный масштаб изменения которой $j_{\mathrm{s}} \gg 1$.

Примечание. Если характереный масштаб изменения $j_{\mathrm{s}}$ последовательности $\left\{a_j\right\}_{j=1}^{\infty}$ составляет хотя бы несколько десятков, ошибку $O\left(1 / j_{\mathrm{s}}^{2 m+2}\right)$ уже при $m=7$ легко уложить в машинную точность.

Метод плавного обрезания#

Нижеописанный метод плавного обрезания немного более ресурсоемок, но зато весьма прост и удобен в программировании, когда необходимо вычислять многомерные бесконечные суммы.

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

Для этого вводится функцию плавного обрезания $j_{\mathrm{sm}}(x)$, которая:

Равна единице при $x<0$.
Равна нулю при $x>1$.
В диапазоне $0 \leq x \leq 1$ она достаточно плавно выходит с единицы на нуль.
Причем касание единицы в $x=0$ и нуля в $x=1$ достаточно плавное.

Приведем два примера такой функции.

Функция

\[\begin{split} j_{\mathrm{sm}}(x)= \begin{cases}\cfrac{1,}{1+\operatorname{expsh} \frac{2 x-1}{X},} & 0 \leq x \leq 1, \\ 0, & x>1,\end{cases} \end{split}\]

в которой величина

\[ X=\frac{1}{\ln \left[\ln \frac{1-\mathrm{e}_{\mathrm{r}}}{\mathrm{e}_{\mathrm{r}}}+\sqrt{1+\ln ^2 \frac{1-\mathrm{e}_{\mathrm{r}}}{\mathrm{e}_{\mathrm{r}}}}\right]} \]

выбирается из условия того, чтобы скачки $j_{\mathrm{sm}}(x)$ в $x=0$ и в $x=1$ заведомо укладывались в машинную точность

\[ j(1)=1-j(0)=\mathrm{e}_{\mathrm{r}} \ll 10^{-16}, \]

причем малая константа $\mathrm{e}_{\mathrm{r}}$ - внутренний параметр метода.

Функция $$ j_{\mathrm{sm}}(x)=\left\{\begin{array}{cc} \left(1-\left(2 x-x^2\right)^{2 k_1}\right)^{k_2}, & 0 \leq x \leq 1 \\ 0, & x>1, \end{array}\right. $$

где $k_{1,2} \in \mathrm{Z}_{+}$, которая обладает касанием высокого порядка единицы при $x=0$ и нуля при $x=1$ и которая вычисляется (при заданных $k_1$ и $k_2$ ) примерно за 8 умножений.

В любом случае функция $j_{\mathrm{sm}}((x-a) /(b-a))$ плавно выходит с единицы при $x<a$ на нуль при $x>b$, где $a<b$.

Опишем теперь сам метод.

Пусть нужно вычислить (многомерную) бесконечную сумму $\sum_{\mathbf{j}} a_{\mathbf{j}}$, где $\mathbf{j}=\left\{j_1, \ldots j_n\right\} \in \mathrm{Z}^n-n$-мерный целочисленный вектор. Считаем, что последовательность $a_{\mathrm{j}}$ при значениях модуля $|\mathrm{j}|$, превышающих некоторое $a \geq 0$, является медленной функцией вектора j: ее характерный масштаб изменения $j_{\mathrm{s}} \gg 1$. Суть метода заключается в следующей формуле

\[ \sum_{\mathbf{j}} a_{\mathbf{j}}=\sum_{\mathbf{j}} a_{\mathbf{j}}\cdot j_{\mathrm{sm}}\left(x_{\mathbf{j}}\right)+\int_{-\infty}^{\infty} d \mathbf{j} \cdot a_{\mathbf{j}} \cdot\left(1-j_{\mathrm{sm}}\left(x_{\mathbf{j}}\right)\right), \]

где функция, например,

\[ x_{\mathbf{j}}=\frac{|\mathbf{j}|-a}{b-a} \text { или } x_{\mathbf{j}}=\frac{|\mathbf{j}|^2-a^2}{b^2-a^2} \text { или др., } \]

а отрезок $|\mathbf{j}| \in[a ; b]$ выбирается так, чтобы хвостовой интеграл выше брался от плавной функции (так что в пределе $j_{\mathrm{s}} \rightarrow \infty$ получается $b-a \sim j_{\mathrm{s}}$.

В этих двух формулах и выполнены следующие трюки:

Последовательность $a_{\mathrm{j}}$ путем домножения на функцию плавного обрезания $j_{\mathrm{sm}}\left(x_{\mathbf{j}}\right)$ с ростом модуля $|\mathbf{j}|$ плавно выводится на нуль, когда $|\mathbf{j}|$ меняется от $a$ до $b$. Так что в многомерной сумме в первом слагаемом в (43) число слагаемых конечно - порядка $b^n$.
В обрезывающем множителе $j_{\mathrm{sm}}\left(x_{\mathrm{j}}\right)$ монотонно возрастающая функция $x_{\mathbf{j}}$ зависит от вектора $\mathbf{j}$ плавно. То же касается последовательности $a_{\mathbf{j}}$ и самого множителя $j_{\mathrm{sm}}\left(x_{\mathbf{j}}\right)$. Это обеспечивает плавность функции под интегралом во втором слагаемом. Поэтому сам интеграл можно численно вычислять схемами старшего порядка, т. е. небольшим компьютерным ресурсом.
Во втором слагаемом в (43) выполнен переход от суммы к интегралу по причине плавности подинтегральной функции. Это законно в силу суммы ряда 2.

Таким образом, суть метода заключается в том, чтобы быстро изменяющуюся составляющую последовательности рассчитывать прямо, а медленно изменяющуюся вычислять путем перехода от суммы к интегралу.

При $m=7$ и $j_{\mathrm{s}} \sim 20$ ошибка $O\left(1 / j_{\mathrm{s}}^{2 m+2}\right)$ перехода от суммы к интегралу в сумме ряда 2, она же - ошибка вычисления интеграла в последней формуле при подобранных $a$ и $b$ - обе вмещаются в машинную точность.

Метод Куммера#

Пусть ряд $\sum_{n=1}^{\infty} a_n$ сходится и сумма его равна S. Предположим, что он сходится крайне медленно и для вычисления его суммы напрямую понадобилось бы слишком много арифметических действий, что приведёт к большой случайной ошибке.

Подберем ряд $\sum_{n=1}^{\infty} b_n$ сумма которого известна $\mathrm{S}_1$, причем, этот ряд в пределе приближает наш ряд т.е. $\lim _{n \rightarrow \infty} \frac{a_n}{b_n}=\gamma \neq 0$

Представим ряд $\sum_{n=1}^{\infty} a_n$ в виде

\[ \sum_{n=1}^{\infty} a_n=\sum_{n=1}^{\infty} a_n-\gamma \cdot \sum_{n=1}^{\infty} b_n+\gamma \cdot \sum_{n=1}^{\infty} b_n=\gamma \cdot S_1+\sum_{n=1}^{\infty} a_n\left(1-\gamma \frac{b_n}{a_n}\right) \]

Так как

\[ \lim _{n \rightarrow \infty} \frac{a_n\left(1-\gamma \cdot \frac{b_n}{a_n}\right)}{a_n}=0 \Leftrightarrow a_n\left(1-\gamma \cdot \frac{b_n}{a_n}\right)=o\left(a_n\right) \]

то ряд $\sum_{n=1}^{\infty} a_n\left(1-\gamma \frac{b_n}{a_n}\right)$ сходится быстрее ряда $\sum_{n=1}^{\infty} a_n$ и поэтому, вычисляя S по формуле

\[ S=\gamma \cdot S_1+\sum_{n=1}^{\infty} a_n\left(1-\gamma \frac{b_n}{a_n}\right) \]

получаем требуемую точность $\varepsilon$ при меньшем числе шагов.

Примеры рядов с известными суммами:

\[\begin{split} \begin{array}{|c|c|} \hline \text { Ряд } & \text { Сумма } \\ \hline \sum_{n=1}^{\infty} \frac{1}{n^2} & \frac{\pi^2}{6} \\ \hline \sum_{n=1}^{\infty} \frac{1}{n^4} & \frac{\pi^4}{90} \\ \hline \sum_{n=1}^{\infty} \frac{1}{n^6} & \frac{\pi^6}{945} \\ \hline \sum_{n=1}^{\infty} \frac{(-1)^{n+1}}{n} & \ln 2 \\ \hline \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n^2} & \frac{\pi^2}{12} \\ \hline \sum_{n=1}^{\infty} \frac{1}{n(n+1) \ldots(n+p)} & \frac{1}{p} \cdot \frac{1}{p !} \\ \hline \sum_{n=1}^{\infty} \frac{\operatorname{cos}(n x)}{n} & -\ln \left(2 \operatorname{sin} \frac{x}{2}\right) \\ \hline \sum_{n=1}^{\infty} \frac{\operatorname{cos}(n x)}{n^2} & \frac{\pi^2}{6}-\frac{\pi x}{2}+\frac{x^2}{4} \\ \hline \sum_{n=1}^{\infty} \frac{\operatorname{sin}(n x)}{n^3} & \frac{\pi^2 x}{6}-\frac{\pi x^2}{4}+\frac{x^3}{12} \\ \hline \sum_{n=1}^{\infty} \frac{\operatorname{sin}(n x)}{n} & \frac{\pi-x}{2} \\ \hline \sum_{n=1}^{\infty}(-1)^{n-1} \cfrac{1}{n} \operatorname{cos}(n x) & \ln \left(2 \operatorname{cos} \frac{x}{2}\right) \\ \hline \sum_{n=1}^{\infty}(-1)^{n-1} \cfrac{1}{n^2} \operatorname{cos}\left(\frac{n \pi x}{3}\right) & \frac{\pi^2}{36}\left(3-x^2\right) \\ \hline \end{array} \end{split}\]

Пример. Вычислить сумму ряда $$ \sum_{n=1}^{\infty} \frac{n^2+1}{n^4+n^2+1} \operatorname{Cos}(2 n) $$

с точностью $\varepsilon=10^{-6}$, применяя метод Куммера.

Определим количество членов, необходимое для обеспечения требуемой точности $\varepsilon=10^{-6}$. $$ \begin{aligned} & \left|R_k\right|=\sum_{n=k+1}^{\infty}\left|\frac{n^2+1}{n^4+n^2+1} \operatorname{Cos}(2 n)\right| \leq \sum_{n=k+1}^{\infty} \frac{n^2+1}{n^4+n^2+1}<\sum_{n=k+1}^{\infty} \frac{n^2+1}{n^2\left(n^2+1\right)}=\sum_{n=k+1}^{\infty} \frac{1}{n^2}< \\ & <\int_{k+1}^{\infty} \frac{d x}{x^2}=\frac{1}{k+1}<10^{-6} \end{aligned} $$

Решая последнее неравенство, получим $$ \mathrm{k}+1>10^{-6}, \quad \mathrm{k}>10^{-6}-1 . $$

Итак, для обеспечения точности $\varepsilon=10^{-6}$ необходимо взять примерно $10^6$ членов ряда - многовато будет.

Теперь воспользуемся вспомогательным рядом из таблицы $$ \sum_{n=1}^{\infty} \frac{\operatorname{cos}(n x)}{n^2}=\frac{\pi^2}{6}-\frac{\pi x}{12}+\frac{x^2}{4} \text { при } \mathrm{x}=2 $$

имеем $$ \sum_{n=1}^{\infty} \frac{\operatorname{cos}(2 n)}{n^2}=\frac{\pi^2}{6}-\frac{\pi}{6}+1 . $$

Находим $$ \gamma=\lim _{n \rightarrow \infty} \frac{a_n}{b_n}=\lim _{n \rightarrow \infty} \frac{\left(n^2+1\right) \operatorname{cos} 2 n \cdot n^2}{\left(n^4+n^2+1\right) \cdot \operatorname{cos} 2 n}=1 $$

Отсюда $$ \begin{aligned} & S=\sum_{n=1}^{\infty} \frac{n^2}{n^4+n^2+1} \cdot \operatorname{cos} 2 n=\frac{\pi^2}{6}-\frac{\pi}{6}+1+\sum_{n=1}^{\infty}\left(\frac{n^2}{n^4+n^2+1}-\frac{1}{n^2}\right) \cdot \operatorname{cos} 2 n= \\ & =\frac{\pi^2}{6}-\frac{\pi}{6}+1-\sum_{n=1}^{\infty} \frac{\operatorname{cos} 2 n}{n^2\left(n^4+n^2+1\right)} . \end{aligned} $$

Определим количество членов полученного ряда, необходимое для обеспечения точности $\varepsilon$. $$ \left|R_k\right|=\left|\sum_{x=k+1}^{\infty} \frac{\operatorname{cos} 2 n}{n^2\left(n^4+n^2+1\right)}\right| \leq \sum_{n=k+1}^{\infty} \frac{1}{n^2\left(n^4+n^2+1\right)}<=\sum_{n=k+1}^{\infty} \frac{1}{n^6}<\int_{k+1}^{\infty} \frac{d x}{x^6}=\frac{1}{5(k+1)^5}<10^{-6} . $$

Откуда находим, что искомая точность достигается при $\mathrm{k} \geq 10$ слагаемых. Можно посчитать и в уме, как говорится.