Рекуррентная формула для синуса

[profrotter]

anik в сообщении #483007 писал(а):

Сравнил, стало интересно. Рекурсия действительно на границе устойчивости, если $k$ отличается от $2\cos\alpha$. Если точнее задать два первых значения, то проблем нет при любой степени дискретизации, причем в диапазоне многих периодов. Моя формула точная.

Не хочу вас расстраивать, но "вашу" формулу получает всякий старательный студент, который изучает курс "дискретная математика", выполняя домашнее задание или, скажем, типовой расчёт, когда решает задачу: Найти общее решение однородного разностного уравнения второго порядка вида:

$$y_n=Ky_{n-1}-y_{n-2}$$

Составляем характеристическое уравнение

$$\lambda^n=K\lambda^{n-1}-\lambda^{n-2}$$

или

$$\lambda^2-K\lambda+1,\lambda\neq 0$$

Корни характеристического уравнения:

$$\lambda_{1,2}=\frac {K\pm\sqrt{K^2-4}} 2$$

В случае, когда $|K|\leqslant 2$, характеристическое уравнение имеет два комплексно-сопряжённых корня:

$$\lambda_{1,2}=\frac {K\pm i\sqrt{4-K^2}} 2$$

Поскольку $|K|\leqslant 2$, можно подобрать такое $\alpha$, что $K=2\cos(\alpha)$. С учётом введённых обозначений:

$$\lambda_{1,2}=\frac {2\cos(\alpha)\pm i\sqrt{4-4\cos^2(\alpha)}} 2=$$

$$=\cos(\alpha)\pm i\sqrt{1-\cos^2(\alpha)}=$$

$$=\cos(\alpha)\pm i\sin(\alpha)=e^{\pm i\alpha}$$

Общее решение исходного разностного уравнения в рассматриваемом частном случае:

$$y_n=C_1\lambda_1^n+C_2\lambda_2^n=C_1e^{in\alpha}+C_2e^{-in\alpha}$$

или (подрбнее этот переход см. topic44178.html)

$$y_n=A\cos(n\alpha)+B\sin(n\alpha),$$

где $C_1,C_2,A,B$ - произвольные постоянные.
Чего теперь хочеться? Синуса? - Тогда полагаем $A=0$. Косинуса? - Тогда $B=0$. А можно получить и рекурсию для более общего случая $y_n=A_m\cos(n\alpha+\varphi)$
Давайте захотим синуса. Неопределённую константу $B$ определяем тогда из начального условия $y_1=\sin(\alpha)$, откуда $B=1$.
Решения рассматриваемого разностного уравнения находится на границе устойчивости не потому, что $K=2\cos(\alpha)$, а потому, что оно является не убывающим на бесконечности и модуль корней характеристического уравнения равен единице. А в случае, если будут иметь место погрешности при вычислении $K$ и оно окажется таким, что $|\lambda_{1,2}|>1$, то вычислительный процесс окажется неустойчивым (это особенно существенно, если вы намерены считать несколько сотен периодов, как вы тут заявляли), другие погрешности вычисления $K$ приведут лишь к тому, что будет вычилсяться рекурсия для угла $\alpha=\arccos(\frac K 2)$ (естественно пока $|K|\leqslant 2$). Если корни характеристического уравнения окажутся не комплексными - то рекурсия даст комбинацию степенных последовательностей.

[INGELRII]

anik в сообщении #482680 писал(а):

Еще могу добавить, формула точная, поэтому погрешность частоты будет зависеть только от точности задания начальных значений синуса и косинуса дискреты угла, и от разрядности вычислителя.

Что формула точная, с этим я спорить не осмеливаюсь. Как бы знаком с такими волшебными словами, как синус суммы, сумма синусов, преобразования эти выполнить способен. Я говорю вот что: у Вас вычисляется рекуррентно. Вы совершаете $n$ шагов. На каждом шаге получается погрешность. Поскольку шагов много, она копится. То есть растет. Весь вопрос в том, насколько быстро! Я дома попытался оценить ее, у меня получилось, что для углов, близких к $0$, она растет аж как степень двойки, то есть экспоненциально. Вот думаю, прав ли я... Но если прав, то формулу можно хоронить, поскольку экспоненциально растущая погрешность - это RIP.

Мои слова, что погрешность становится близка к единице, Вы тоже весьма творчески интерпретировали. Периоды изменения аргумента тут ни при чем. Погрешность же не в значении периода, а в значении собственно самой функции от этого аргумента. Ну вот вычислили Вы, допустим, синус угла $\frac{\pi}{4}$ с погрешностью $1$, и получили, что он равен $1,7$ . И что теперь, учитывая, что синус вещественного числа по модулю не превосходит единицы?

Это ведь проверить легко. Так и быть, займусь дома, хотя лень невероятно. Вычислю по Вашей рекуррентной формуле синус от угла $\frac{\pi}{2}$, разбив его на 10000 частей. Посмотрим, с какой погрешностью он мне его выдаст. Но предчувствия у меня самые нехорошие.

[anik]

Уважаемый profrotter! Спасибо за информацию. Я познакомлюсь с разностными уравнениями, потом Вам отвечу.

[anik]

INGELRII в сообщении #483245 писал(а):

Вычислю по Вашей рекуррентной формуле синус от угла $\frac{\pi}{2}$, разбив его на 10000 частей. Посмотрим, с какой погрешностью он мне его выдаст. Но предчувствия у меня самые нехорошие.

У меня тоже нехорошие предчувствия.
Если $\frac{\pi}{2}$ поделить на 10000 частей, получится угол около 30 угловых секунд. Необходимо взять 14 знаков после запятой в синусе и косинусе такой дискреты угла, чтобы получить четыре значащих цифры. А вообще в чём практическая необходимость брать такие маленькие дискреты угла, чтобы доказать, что точность поплывёт?

[INGELRII]

Ну у меня есть стойкое убеждение, что она поплывет. Я хочу его обосновать. Чтобы сделать это наверняка, я и беру такие запредельно малые дискреты... Так сказать overkill. Кстати, если на таких малых углах точность не поплывет, это будет уже серьезный аргумент в Вашу пользу. Но отчет об этом будет завтра, я с работы пишу. Все математические программы у меня дома. Терпение.

[profrotter]

 int N = 40000;
 double x1 = 0.0, x2 = 2000.0*M_PI;
 double a = x2 / (double)(N-1);
 float y[42000];
 float K = 2.0*cos(a);
 y[0]=sin(0.0);
 y[1]=sin(a);
 for(int n = 2; n < N; n++)
  y[n] = K*y[n-1]-y[n-2];

40000 итераций, тысяча периодов синуса при шаге дискретизации $\alpha= \frac {\pi} {20}$. float - точности хватает.

Сплошной - синус из Borland C++.

[anik]

profrotter, спасибо за Ваш труд.

[arseniiv]

По графику не видно некоторых вещей. Вы не могли бы сделать отдельно график разностей?

[profrotter]

arseniiv в сообщении #483551 писал(а):

По графику не видно некоторых вещей. Вы не могли бы сделать отдельно график разностей?

Теперь уже до понедельника - программка на работе осталась. Удивительно вот что: допустим на каком-то шаге ввиду конечной разрядности сетки отсчёт оказался вычисленным с некторой точностью, то есть $y_n=\bar{y}_n+\Delta$. Тогда на следующем шаге эта погрешность превращается в $K\Delta$, потом в $K^2\Delta$ и т.д. Это означает, что при $|K|>1$ или $\cos(\alpha)>\frac 1 2$ или $\alpha < 60^{\circ}$ будет иметь место накопление погрешности, причём чем меньше шаг дискретизации, тем ближе $\cos(\alpha)$ к единице, а значение $|K|$ к 2. Только не сама погрешность нарастает экспоненциально, как писал INGELRII, а множитель, на который умножается эта погрешность на каждом последующем шаге. Несмотря на это существенных расхождений я не получил.
Даже если максимальное отклонение окажется порядка, скажем $10^{-3}$, формулу можно будет признать пригодной для генерации гармоник, например, в том случае, когда результаты будут восприниматься на слух, на что и указвал выше автор темы. Интересно было бы также посмотреть что получится, если ограничить точность вычислений, например, 4-мя знаками после запятой. Доживём до понедельника. А вообще-то задачка простая и форум большой и только profrotter сидит и строит графики...

[arseniiv]

(Оффтоп)

profrotter в сообщении #483574 писал(а):

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

Неужели он у людей такой плохой?

[anik]

arseniiv в сообщении #483598 писал(а):

(Оффтоп)

profrotter в сообщении #483574 писал(а):

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

Неужели он у людей такой плохой?

При настройке роялей пользуются камертоном, частота колебания которого 440 Гц (Ля, первой октавы). Струну настраивают по частоте биений. Если период биений 4-5 сек, то такой точности достаточно. Это уход частоты на один период примерно за 2000 колебаний. Для гармоник звукового спектра тембра, такая точность, наверное, не требуется.

[ewert]

anik в сообщении #482941 писал(а):

Эта формула мне известна уже давно.

Не только Вам. В алгоритме быстрого преобразования Фурье испокон веков именно так и поступают (правда, в комплексном варианте, но какая разница).

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

Ну и по применению. Это стандартный приём при вычислении БПФ, в алгоритме Герцля (определяющим, есть ли в сигнале определённые частоты, например, при дешифровке телефонного тонового набора) и т.д. где нужны последовательные значения синуса (и косинуса, их обычно считают одновременно)
Проще всего вывести из тождества $e^{ix}=\cos x + i \sin x$

[anik]

Евгений Машеров в сообщении #483712 писал(а):

Проще всего вывести из тождества $e^{ix}=\cos x + i \sin x$

Интересно было бы посмотреть на вывод рекуррентной формулы для синуса из формулы Эйлера, тем более, что это "проще всего вывести".

[ewert]

anik в сообщении #483726 писал(а):

Интересно было бы посмотреть на вывод рекуррентной формулы для синуса из формулы Эйлера, тем более, что это "проще всего вывести".

$e^{i(k+1)h}=e^{ikh}e^{ih}=(\cos kh+i\sin kh)(\cos h+i\sin h)$ и т.д.

Это не то что бы проще, но не нужно помнить тригонометрию.