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

[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 (остальное отрезал)

[ihq.pl]

wrest в сообщении #1726979 писал(а):

Самое время попробовать произвольную точность в «интерактивный курс: введение в программирование на PARI/GP»

Да, давно пора, спасибо
Вы в PARI/GP запрограммировали эту формулу? Наверняка, там можно сделать это проще и нагляднее.

wrest в сообщении #1726979 писал(а):

Сумма будет 0.8813507609

У меня при тех же значениях получилось 0.867137. Причина такого сильного расхождения явно не в точности(
И снова похожая картина: значение 0.867137 достигается при n=12 и больше не меняется до n=100 (на этом уровне точности). А мой первый график со скачками - это то, что получилось при

$$\frac{x_0^2}{r^2}=1/3,\quad \frac{y_0^2}{r^2}=1/6,\quad \frac{a}{r}=0.9\cos(\pi/3),\quad \frac{b}{r}=0.9\sin(\pi/3),\quad \pi=3.14159265358979.$$

[wrest]

ihq.pl
Вот код:

Sn_vector(a, b, r, x0, y0, n) = {
    my( a2 = a^2 / r^2,
        b2 = b^2 / r^2,
        x2 = 4 * x0^2 / r^2,
        y2 = 4 * y0^2 / r^2,
        maxFact = 4*n,
        F = vector(maxFact + 1, i, (i-1)!),
        coeff = 4 * a * b / r^2,
        
        Omega_cache = matrix(n+1, n+1),
        
        current_S,
        S = vector(n+1, i, 0),
        delta );

    Omega_cache[1, 1] = 1;
    current_S = coeff;
    S[1] = current_S;

    for(n = 1, n,
        delta = 0;
        for(m = 0, n,
            my( n_idx = n - m );
            
            \\ Вычисляем Omega сразу, без проверки
            Omega_cache[m+1, n_idx+1] = sum(j = 0, m,
                (-1)^j * x2^j *
                sum(i = 0, n_idx,
                    (-1)^i * y2^i *
                    F[m+j+n_idx+i+1] / (F[m-j+1] * F[n_idx-i+1] * F[2*j+1] * F[2*i+1])
                )
            );
            
            delta += (-1)^n / ((2*m + 1) * (2*n_idx + 1)) *
                     a2^m * b2^n_idx * Omega_cache[m+1, n_idx+1]
        );
        current_S += coeff * delta;
        S[n+1] = current_S
    );

    return(S);
}

Вставляете его в синий бокс тут: https://pari.math.u-bordeaux.fr/gp.html
Нажимаете "evaluate with pari", и затем нажимаете "clear input"
Теперь копируете строчку
v=Sn_vector(1,1,2,1/2,1/3,50); print("Ready");
в тот же синий бокс и нажимаете "evaluate with pari"
Ждёте
В красном боксе появится Ready
Теперь в синем боксе печатаете
v[45]+.
Нажимаете "evaluate with pari"
В красном боксе появлется Sn для n=45
Ну и т.п.
Аргументы: Sn_vector(a, b, r, x0, y0, n)
Функция возвращает вектор из n рациональных чисел (точных!) соответствующих частичным суммам вашего ряда.
Еслииаргументы рациональные, и если напечатать в синий бокс v[5] (то есть без +.) то получите точное значение в виде p/q
? v[5]
%6 = 596763551383/677221171200

Факториалы и $\Omega_{m,k}$ кешируются чтобы много раз не вычислять, ну и за один проход вычисляются все частичные суммы, по скорости практически то же самое что вычислить последнюю.