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

[anik]

ewert в сообщении #483742 писал(а):

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

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

Извините, но я ничего не понял. Наверное, туповат.

[ewert]

Д.т.:

$\cos(k+1)h+i\sin(k+1)h=\cos kh\cdot\cos h+i\sin kh\cdot\cos h+i\sin h\cdot\cos kh-\sin kh\cdot\sin h;$

$\begin{cases}\cos(k+1)h=c\cdot\cos kh-s\cdot\sin kh,\\ \sin(k+1)h=s\cdot\cos kh+\cdot\sin kh,\end{cases}$

(где, естественно, $c\equiv\cos h$ и $s\equiv\sin h$ вычисляются раз и навсегда, а дальше уже только арифметика). Только видите ли в чём дело: если речь о спектральном анализе, то программировать в синусах и косинусах просто неудобно -- гораздо естественнее реализовать всё это в комплексной арифметике.

[INGELRII]

profrotter
Уже смотрел. Да, растет не сама погрешность, а множитель.

all
И все равно. Считал по формуле $\sin(\frac{\pi}{2})=1$, количество итераций $10000$, начальная точность $10$ знаков после запятой. На последнем получилось $0,9934535098$, что от единицы далековато. Причем на шаге с номером $9937$ "синус" начал уже убывать, а максимальное значение так и не достигло единицы. Че-то как-то...

Не знаю, может кому и нравятся формулы, которые десять верных знаков превращают в два после ${10}^{4}$ шагов. Но я бы взял в руки лопату и выкопал могилку для формулки.

[anik]

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

Не знаю, может кому и нравятся формулы, которые десять верных знаков превращают в два после ${10}^{4}$ шагов.

Я думал, что точность выражается не в количестве десятичных знаков после запятой, а в количестве значащих цифр. Может быть я ошибался?
Если за дискрету угла взять секунду, то шесть знаков после запятой для синуса будут нули, а для косинуса - 10 девяток. Какова при этом будет точность?

[profrotter]

Смотрим на разность.
По всем вычислениям:

Последние пара периодов:

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

-- Пн сен 19, 2011 17:43:43 --

Кстати, вот и сама программка http://depositfiles.com/files/6zrf4skgd

[anik]

Правильно, я так и поступал. Обновлял вычисления, когда датчик угла обнулялся.

[profrotter]

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

И все равно. Считал по формуле $\sin(\frac{\pi}{2})=1$, количество итераций $10000$, начальная точность $10$ знаков после запятой. На последнем получилось $0,9934535098$, что от единицы далековато. Причем на шаге с номером $9937$ "синус" начал уже убывать, а максимальное значение так и не достигло единицы. Че-то как-то...

Позвольте, но при угле $\frac {\pi} 2$ параметр $K=2\cos(\frac {\pi} 2)=0$ и рассматриваемая рекуррентная формула принимает вид: $y_n=-y_{n-2}$. Как вы умудрились получить в этом случае погрешность?!

[INGELRII]

profrotter
Да легко. Как следует из моего сообщения, число $\frac{\pi}{2}$ я разделил на ${10}^{4}$. Получил угол $\frac{\pi}{20000}$ (косинус от него не ноль, а довольно близок к единице), который и подставил в горячо обсуждаемую нами формулу. Сделал по ней ${10}^{4}$ итераций, и на последней получил "синус" от того числа, которое получается умножением $\frac{\pi}{20000}$ на ${10}^{4}$. Ошибка, повторюсь, получается уже в третьем знаке. Буду рад дать еще более подробные пояснения, коли в них возникнет нужда.

Как Вы понимаете, при этом я вообще ни на какие тысячные периоды не выходил. Вообще не сделал ни одного периода. Я путешествовал от нуля до пи-пополам. Синус вычислил от пи-пополам.

Тестовые задачи меня учили выбирать так, чтобы они удовлетворяли простому критерию. Не нужно смотреть, сколько там получается синус после миллиона периодов, ни к чему вся эта мегаломания. Если мы тестируем формулу, то для начала надо взять какой-нибудь примитивный случай. Скажем, все знают, что $\sin(\frac{\pi}{2})=1$. Вот и подберем значение угла и число итераций так, чтобы на финальном шаге получался синус от пи-пополам! Так я и сделал. И только в том случае, если формула показала себя достойно, мы начинаем по ней считать $\sin({10}^{800}\pi)$ ради прикола. А не в обратном порядке.

anik
Я брал первые десять значащих цифр. То есть ненулевых. Согласен насчет косинуса. там получается слишком маловато не-девяток в конце. Хорошо, пересчитаю заново. Возьму точность не десять, а двадцать значащих цифр. Этого будет достаточно?

[profrotter]

INGELRII, не могли бы вы привести программный код?

-- Вт сен 20, 2011 13:03:37 --

У меня в описанном вами случае получается совсем другое:

Всё выглядит так, как будто бы вычислительная процедура вообще теряет устойчивость. Это и не удивительно: ссылку на отдельную тему, где рассматриваются вопросы вычисления синуса я уже в этой теме давал, там же и указал на взможность применения данной формулы для расчётов, проблемах устойчивости при малых шагах дискретизации, написал о ограничениях и достигаемой точности.
Собственно, не понятно почему вы спешите хоронить форумулу, а не хотите задуматься о границах её применения? Почему не учитываете сообщения ewert и Евгения Машерова? Речь идёт о формировании сигналов при цифровой обработке. Там важнее скорость вычислений, а точность порядка $10^{-3}$ является вполне приемлимой (конечно в зависимости от конкретно решаемой задачи) и уж конечно является нелепым синтезировать гармонический сигнал, дискретизируя его с шагом $\frac {\pi} {20000}$. Гармонический сигнал можно восстановить, когда шаг дискретизации не более $\frac {\pi} 2$. На практике шаг дискретизации должен выбираться в несколько (5-10) раз меньше.

-- Вт сен 20, 2011 13:20:57 --

Смотрим расчёт синуса при шаге дискретизации $0,1$ с точностью double:

То же с округлением до 4 знаков:

[profrotter]

К сожалению не смог отредактировать - в послдеднем сообщение картинка для точности float.

[anik]

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

Возьму точность не десять, а двадцать значащих цифр. Этого будет достаточно?

Если вы возьмёте двадцать десятичных знаков после запятой, то для косинуса это будет 13 значащих цифр, а для синуса 17 значащих цифр. (Здесь речь идёт о начальных значениях функций дискреты угла). Этого более чем достаточно для большинства практических потребностей.

[profrotter]

Тема пережила этап своего активного обсуждения и скоро найдёт себе достойное место в архивах форума. Пока этого окончательно не произошло я хотел немного рассказать об одном методе построения рекуррентных формул для различных функций.
Искомую линейную рекуррентную формулу запишем в виде:

$$y_n=\sum\limits_{k=0}^{K}a_k x_{n-k}+\sum\limits_{m=1}^{M}b_m y_{n-m}\eqno (1)$$

Построение формулы заключается в подборе таких коэффициентов $\{a_n\}$ и $\{b_m\}$ и такой дискретной последовательности $x_n$, что результат расчёта будет давать отсчёты некоторой заданной функции $y_n=f(nT)$, $T$ - период дискретизации.
В общем случае такая задача может иметь большое количество решений, а само нахождение их плохо поддаётся формализации. Рассмотрим частный случай, когда

$$x_n=\delta_n=\begin{cases}1,&n=0\\0,&n\neq 0\end{cases} \eqno (2)$$

Заметим, что в рассматриваемом случае следующие равенства дают начальные условия для искомой рекуррентной формулы:

$$y_0=a_0$$

$$y_1=a_1+b_1 y_0$$

$$y_2=a_2+b_1 y_1+b_2 y_0$$

$$y_3=a_3+b_1 y_2+b_2 y_1+b_3 y_0$$

$$\cdots$$

$$y_K=a_K+b_1 y_{K-1}+b_2 y_{K-2}+\cdots +b_K y_0\eqno (3)$$

А вместо (1) можно записать

$$y_n=\sum\limits_{m=1}^{M}b_m y_{n-m}, &n>K\eqno(4)$$

Вернёмся к (1) и возмём Z - преобразование от обеих частей равенства и, с учётом свойства временного запаздывания и того, что $x_n\leftrightarrow 1$, получим:

$$Y(z)=\frac {\sum\limits_{k=0}^{K}a_k z^{-k}} {1-\sum\limits_{m=1}^{M}b_m z^{-m}}\eqno (5)$$

Теперь пострение рекуррентного соотношения можно предствить в виде следующей последовательности этапов:
1. Осуществляем дискретизацию $f(t)$ и получаем последовательность $y_n=f(nT)$;
2. Находим Z - преобразование полученной последовательности и приводим к виду (5);
3. Находим коэффициенты $\{a_n\}$ и $\{b_m\}$;
4. Подставляем коэффициенты либо в (1) с последовательностью (2), либо в (4) с начальными условиями (3).
Смотрим примеры:
1. Экспоненциальный импульс $f(t)=\sigma(t)e^{-\alpha t}$, дискретизируем $y_n=\sigma(nT)e^{-\alpha nT}$. В таблице Z - преобразований находим $Y(z)=\frac 1 {1-e^{-\alpha T}z^{-1}}$, откуда $a_0=1,b_1=e^{-\alpha T}$. Рекуррентные формулы:

$$y_n=\delta_n+e^{-\alpha T} y_{n-1}$$

или

$$y_n=e^{-\alpha T} y_{n-1}, &n>0$$

с начальным условием $y_0=1$.
2. Экспоненциальный радиоимпульс

$$f(t)=\sigma(t) e^{-\alpha t} \sin(\omega_0 t)$$

$$y_n=\sigma(nT) e^{-\alpha nT} \sin(\omega_0 nT)$$

$$Y(z)=\frac {e^{-\alpha T}\sin(\omega_0 T)z^{-1}} {1-2e^{-\alpha T}\cos(\omega_0 T)z^{-1}+e^{-2\alpha T}z^{-2}}$$

$$a_0=0,a_1=e^{-\alpha T}\sin(\omega_0 T),b_1=2e^{-\alpha T}\cos(\omega_0 T),b_2=-e^{-2\alpha T}$$

Рекуррентные формулы:

$$y_n=e^{-\alpha T}\sin(\omega_0 T)\delta_{n-1}+2e^{-\alpha T}\cos(\omega_0 T)y_{n-1}-e^{-2\alpha T}y_{n-2}$$

или

$$y_n=2e^{-\alpha T}\cos(\omega_0 T)y_{n-1}-e^{-2\alpha T}y_{n-2},&n>1$$

с начальными условиями

$$y_0=0,y_1=e^{-\alpha T}\sin(\omega_0 T)$$

В частном случае при $\alpha=0$ получим формулы для синуса:

$$y_n=\sin(\omega_0 T)\delta_{n-1}+2\cos(\omega_0 T)y_{n-1}-y_{n-2}$$

или

$$y_n=2\cos(\omega_0 T)y_{n-1}-y_{n-2},&n>1$$

с начальными условиями

$$y_0=0,y_1=\sin(\omega_0 T)$$

[anik]

profrotter, ну Вы молодец!
Я эту формулу нашёл выпендриваясь с микрокалькулятором, когда они ещё только начали появляться в СССР, чисто случайно. Вы вывели её из общей теории. Этот метод достоин того,чтобы его изучить.
Ещё раз, большое Вам спасибо!

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

Это было бы ещё замечательнее, если бы не было бы в любом учебнике по цифровым фильтрам...

[profrotter]

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

profrotter, ну Вы молодец!
Я эту формулу нашёл выпендриваясь с микрокалькулятором, когда они ещё только начали появляться в СССР, чисто случайно. Вы вывели её из общей теории. Этот метод достоин того,чтобы его изучить.
Ещё раз, большое Вам спасибо!

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

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

Это было бы ещё замечательнее, если бы не было бы в любом учебнике по цифровым фильтрам...

Совершенно с вами согласен: именно от цифровых фильтров я и отталкивался. По сути это задача синтеза рекурсивного цифрового фильтра по аналоговому фильтру-прототипу методом инвариантной импульсной характеристики. Затем подаём на вход цифрового фильтра единичный отсчёт и получаем что требовалось. Однако, во-первых в теме никто не упомянул о связи рекуррентных формул и теории линейных дискретных систем (я показал автору вывод формулы путём решения уравнения, которое свалилось будто бы с потолка, вы вместе с ewert дали доказательство через комплексные числа - в обоих случаях формализованный путь к получению формул не был указан), во-вторых в учебниках есть материал по разностным уравнениям, Z-преобразованию, цифровым фильтрам, но нигде не формализована последовательность этапов получения именно форумул для расчёта значений функций (если я тут не прав, то прошу подсказать мне учебник, где излагается "рецепт", похожий на мой ), в-третьих я ориентировался на то, что большинство участников темы (в частности автор, который признался, что не знаком с разностными уравнениями) не знакомы с цифровыми фильтрами, потому пытался изложить материал не прибегая к понятию цифрового фильтра, импульсной характеристики и тп. Собственно, я и не претендовал на какую-либо уникальность и новизну - так и написал, что "хочу рассказать...".