Сходимость степенного ряда

[ihq.pl]

У меня есть функция

$$g(u,v) = [(x_0+u)^2+(y_0+v)^2 + c^2]^{-1}.$$

Ряд Маклорена этой функции сходится в круге

$$u^2+v^2<r^2, \quad r^2= x_0^2 + y_0^2 + c^2.$$

Поэтому

$$\int_{-a}^a\int_{-b}^b g(u,v)dudv = \sum_{n=0}^\infty \sum_{m=0}^n g_{n,m}\int_{-a}^a\int_{-b}^b u^mv^{n-m}dudv,\quad a^2+b^2<r^2.$$

Выбираю подходящие $a,b,x_0, y_0$ и считаю сумму первых $n$ членов полученного ряда. Сначала всё хорошо: начиная с некоторого малого $n$ эта сумма почти не меняется. Но начиная с некоторого большого $N$ сумма начинает резко расти. Причем, буквально для $n=N+10$ эта сумма принимает запредельные значения. Перепроверял сто раз. Вроде ошибок нет. Может, это как-то связано с так называемой нестабильностью?

[ihq.pl]

ihq.pl в сообщении #1726839 писал(а):

$$\int_{-a}^a\int_{-b}^b g(u,v)dudv = \sum_{n=0}^\infty \sum_{m=0}^n g_{n,m}\int_{-a}^a\int_{-b}^b u^mv^{n-m}dudv,\quad a^2+b^2<r^2.$$

Лучше, наверное, так написать

$$\int_{-a}^a\int_{-b}^b g(u,v)dudv = S_n + R_n$$

где

$$S_n = \sum_{k=0}^n \sum_{m=0}^k g_{k,m}\int_{-a}^a\int_{-b}^b u^mv^{k-m}dudv, \quad g_{k,m} = \frac{1}{k!}C_k^m\frac{\partial^k g}{\partial u^m\partial v^{k-m}}\Big|_{(0,0)}.$$

Последовательность $S_n$ должна по идее сходиться, но ведет себя странно. Быстро приближается к некоторому числу. Для $6\leqslant n\leqslant 45$ почти не меняется, и потом начинает резко расти.

[mihaild]

А Вы это всё проверяете в арифметике произвольной точности, или с плавающей запятой?

[Евгений Машеров]

Это скорее всего чисто вычислительная ошибка. Я бы, исходя из "стремительного роста", поискал бы ошибку в расчёте коэффициентов, что-то с факториалами.

[ihq.pl]

mihaild, с плавающей запятой (8 байт).

-- добавлено через 4 минуты --

Евгений Машеров, в точности кодирую в С эту формулу:

$$S_n = \\ &\frac{4ab}{r^2}\sum_{m=0}^{n}\frac{(-1)^m}{2m+1}\Big(\frac{a^2}{r^2}\Big)^m \sum_{k=0}^{n-m}\frac{(-1)^k}{2k+1}\Big(\frac{b^2}{r^2}\Big)^k\Omega_{m,k},$$

$$\Omega_{m,k} = \sum_{j=0}^m (-1)^j 4^j\Big(\frac{x_0^2}{r^2}\Big)^j \sum_{i=0}^k (-1)^{i}4^i \Big(\frac{y_0^2}{r^2}\Big)^iC_{m,k}^{i,j},\quad C_{m,k}^{i,j} = \frac{(m+j+k+i)!}{(m-j)!(k-i)!(2j)!(2i)!}.$$

Это - явная формула для $S_n$ из моего предыдущего сообщения. Cделал только совсем маленькое преобразование:

$$\sum_{k=0}^n\sum_{m=0}^k g_{k,m}\int u^mv^{k-m}dudv = \sum_{m=0}^n \sum_{k=0}^{n-m} g'_{k,m}\int u^mv^kdudv, \quad g'_{k,m} = \frac{1}{(k+m)!} C_{k+m}^m \frac{\partial^{k+m}}{\partial u^m\partial v^k}\Big|_{(0,0)}.$$

Самым сложным здесь было для меня вычисление коэффициентов $C_{m,k}^{i,j}$. Проверял вручную для $n=2$. Вычисляются итеративно. Не так чтоб сначала один факториал, потом другой...

-- добавлено через 47 минут --

График $S_n$. Дальше стремительный рост по амплитуде.

[ihq.pl]

ihq.pl в сообщении #1726940 писал(а):

Вычисляются итеративно

Рекурсивно: в циклах по $i$ и $j$ результат предыдущего шага умножается на соответствующий множитель.

[Someone]

ihq.pl в сообщении #1726940 писал(а):

с плавающей запятой (8 байт).

ihq.pl в сообщении #1726948 писал(а):

в циклах по $i$ и $j$ результат предыдущего шага умножается на соответствующий множитель.

Почти наверняка Вы наблюдаете результат накопления ошибок округления.

Попробуйте сделать вычисления с двойной точностью. Эффект будет тот же, но при бо́льших значениях $n$.
Это хорошо известное явление.

[ihq.pl]

Someone в сообщении #1726957 писал(а):

Почти наверняка Вы наблюдаете результат накопления ошибок округления.

Попробуйте сделать вычисления с двойной точностью. Эффект будет тот же, но при бо́льших значениях $n$.
Это хорошо известное явление.

К сожалению, мой компилятор не увеличивает точность, когда меняешь тип double на long double. Но я уменьшил. На графике то, что получилось с типом float:


Спасибо.

[dsge]

$r$, $a$, $b$ могут быть больше единицы?
Посмотрите значения $g_{k,m}$ при $k>45$,

[ihq.pl]

dsge в сообщении #1726969 писал(а):

$r$, $a$, $b$ могут быть больше единицы?

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

$$a^2+b^2<r^2,\quad r^2 = x_0^2+y_0^2+c^2.$$

Дроби $x_0/r, y_0/r$ всегда строго меньше единицы т.к. $c>0$. Однако при

$$a=b=1, x_0=y_0=c=1$$

ничего плохого не происходит даже при $n=100$. Вот график при этих значениях:



Значение $S_n = 0.980587$ достигается при $n=8$ и дальше вообще не меняется вплоть до $n=100$. А дальше я не проверял. А почему Вы спросили про значения именно этих параметров?

-- добавлено через 14 минут --

То же самое при

$$a=b=2, x_0=y_0=c=2.$$

[dsge]

Ряд Тейлора может плохо себя вести, например для функции $e^{100x}$. Но этот ряд плохо себя ведет на начальных членах, а у вас как-то странно, начиная с 45-го. Скорее всего в самом деле с этого места факториалы начинают убивать точность, и двойная точность не поможет против $45!$

[ihq.pl]

dsge в сообщении #1726976 писал(а):

с этого места факториалы начинают убивать точность

Это то, что имел в виду Евгений Машеров, но я не понял.

dsge в сообщении #1726976 писал(а):

и двойная точность не поможет против $45!$

К сожалению, не могу пока это проверить.

[wrest]

ihq.pl
Самое время попробовать произвольную точность в «интерактивный курс: введение в программирование на PARI/GP»
Для рациональных чисел будет вообще без потери точности (т.е. точно). Но выглядеть будет страшненько. Но это если надо, конечно. Можно и в приближенном виде с нужным количеством знаков печатать.

для

a  = 1;
b  = 1;
r  = 2;
x0 = 1/2;
y0 = 1/3;
n  = 100;

Сумма будет 0.8813507609 (остальное отрезал)