Рекурсивное получение синусоиды

[Lesobrod]

Необходимо сделать несложную рекурсию, генерирующую синусоиду. Начал с родного уравнения
$ \frac {d^2x}{dt^2}+\omega^2 x=0$
Используя простейшую аппроксимацию производной,
получим рекурсию
$x(k)=\frac {2 x(k-1)-x(k-2)} {1+\omega^2 \Delta t^2},x(0)=x_0,x(1)=x_1 $
Численный расчет показывает, что она даёт затухающие колебания.
Затухание, конечно, тем меньше, чем меньше шаг по времени,
т.к. это явно ошибка аппроксимации.
Я знаю,что существуют гораздо более продвинутые методы решения дифуров,
всякие "коррекции-предсказания" и т.п. Но хочется чего-то простого и чистого....
Как например, вот это:
$ \left\{ \begin{array}{l} x(k)=c x(k-1) + s y(k-1),\\ y(k)=-s x(k-1)+c y(k-1), \end{array} \right $
где $ с=cos(\omega \Delta t), s=sin(\omega \Delta t) $
Такая простая рекурсия даёт отличную незатухающую синусоиду.
Эта система попалась мне в одной статье по цифровому синтезу звука;
откуда её взяли - неизвестно.
С другой стороны, если я правильно пониманию, систему двух разностных уравнений
первого порядка можно свести к одному уравнению второго порядка.
Но что-то не очень получается..
Может, кто-нибудь сможет внести ясность?
Хотелось бы не просто готовой формулы, а неких пояснений

[ewert]

Lesobrod в сообщении #364335 писал(а):

Численный расчет показывает, что она даёт затухающие колебания.

Естественно -- поскольку Ваша двухшаговая схема несимметрична:

$\dfrac{x_{k+1}-2x_{k}+x_{k-1}}{\tau^2}+\omega^2x_{k+1}=0\,.$

Корни соответствующего характеристического уравнения: $(1+\tau^2\omega^2)\lambda^2-2\lambda+1$ по модулю не равны единице (чуть меньше единицы) -- из-за того, что не совпадают первый и последний коэффициенты, отсюда и затухание. С уменьшением временного шага $\tau$ затухание уменьшается, т.к. декремент затухания зависит от $\tau$ квадратично, в то время как частота -- линейно.

А вот если выбрать симметричную схему ("с перешагиванием"):

$\dfrac{x_{k+1}-2x_{k}+x_{k-1}}{\tau^2}+\omega^2x_{k}=0\,$

-- то всё будет нормально. Теперь $\lambda^2+(-2+\tau^2\omega^2)\lambda+1=0$, откуда $\lambda_{1,2}=e^{\pm i\varphi}$, где $\cos\varphi=1-{1\over2}\tau^2\omega^2$ и, значит, при малых шагах $\varphi\approx\tau\omega$ (но всё-таки не совпадает с ним). Чтобы получить точную частоту, надо соотв. образом подкорректировать параметр $\omega$ в уравнении.

Lesobrod в сообщении #364335 писал(а):

Такая простая рекурсия даёт отличную незатухающую синусоиду.Эта система попалась мне в одной статье по цифровому синтезу звука;откуда её взяли - неизвестно.

Эта схема в матричной записи имеет вид:

$\begin{pmatrix}x\\y\end{pmatrix}_{k}=\begin{pmatrix}c&s\\-s&c\end{pmatrix}\cdot\begin{pmatrix}x\\y\end{pmatrix}_{k-1}$

Собственные числа такой матрицы -- это $e^{\pm i\varphi}$, где $\cos\varphi=c$ и $\sin\varphi=s$, что даёт незатухающие колебания с частотой $\dfrac{\varphi}{\tau}$, откуда и те соотношения на $c$ и $s$.

[Lesobrod]

Спасибо, но всё таки из каких физических соображений получена
эта система и как корректно привести её к одному уравнению?

[ewert]

Lesobrod в сообщении #364344 писал(а):

Спасибо, но всё таки из каких физических соображений получена
эта система и как корректно привести её к одному уравнению?

Дались Вам эти физические соображения. Ну запишите, например, исходный дифур второго порядка стандартно как систему двух дифуров первого порядка. Выпишите для неё две одношаговых разностных схемы -- явную и неявную. И скомбинируйте их так, чтобы матрица в правой части имела единичный детерминант. (Скорее всего, их надо просто усреднить с последующим преобразованием к стандартному виду, т.к. это делает схему симметричной.)

Что касается сведения к одному уравнению. Прежде всего -- этого делать не следует. А если уж так хочется, то можно сделать ещё один шаг:

$\begin{pmatrix}x\\y\end{pmatrix}_{k+1}=\begin{pmatrix}c_2&s_2\\-s_2&c_2\end{pmatrix}\cdot\begin{pmatrix}x\\y\end{pmatrix}_{k-1}$

(последняя матрица получается возведением в квадрат первой; естественно, $c_2=c^2-s^2$ и $s_2=2sc$). Теперь выпишите первые строчки этого и исходного векторных равенств:

$\begin{cases}x_{k+1}=c_2x_{k-1}+s_2y_{k-1};\\x_{k}=cx_{k-1}+sy_{k-1}.\end{cases}$

И исключите отсюда $y_{k-1}$. После всех сокращений получится $x_{k+1}-2c\,x_{k}+x_{k-1}=0$ (естественно; а что ещё и получиться-то могло бы).

[Lesobrod]

Всё, ура!
Последний вопрос (если заглянете) ;)))
Хочу ввести нелинейность и взаимодействие нескольких
осцилляторов. Правильно ли я понял, что симметричная схема,
ведущая себя корректно, должна выглядеть примерно так
(для двух осцилляторов):

$ \left\{ \begin{array}{l} x_{k+1}+x_{k-1}+f(x_k)+g(y_k)=0,\\ y_{k+1}+y_{k-1}+f(y_k)+g(x_k)=0, \end{array} \right. $
,
где $f$отвечает за нелинейность, а $g$-за взаимодействие ?

[ewert]

Ну что-то типа.

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

Ошибка тут возникает оттого, что Вы записали модель в виде дифференциальных уравнений, а моделируете конечно-разностными. Которые могут рассматриваться, как аппроксимация дифференциальных при разностях, стремящихся к нулю, но всё же не тождественны им.

[ewert]

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

Ошибка тут возникает оттого, что Вы записали модель в виде дифференциальных уравнений, а моделируете конечно-разностными.

Не из-за этого. У него ведь проблемы были не с неточностью решений, а с нарушением периодичности. А это объясняется не конечноразностностью, а неудачным выбором модели.

Бывают и вовсе удивительные случаи. Скажем, для простейшей задачи Штурма-Лиувилля и её простейшей конечноразностной аппроксимации собственные функции в точности совпадают с истинными (собственные числа, разумеется, лишь приближённо).