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

[Lesobrod]

Для синтеза звука необходимы (в каждый момент времени) значения функции
$f(t) = t^{a}\cdot e^{-\delta t}$
Вычислять их прямо каждую $1/44100$ долю секунды слишком жирно.
Вот решил составить рекуррентную формулу, чтобы по 1-2 предыдущим
значениям получать очередное.
Взяв производную по времени, видим, что:
$df(t)/dt = (a/t - \delta)\cdot f(t)$
Заменяя производную разностным соотношением и учитывая, что $t_n = n dt$, получим:
$(f_n - f_{n-1})/dt = (a/t_n - \delta)\cdot f_{n-1}$
$f_n = (1 + a/n - \delta\cdot dt)\cdot f_{n-1}$
Такой процесс работает, но при больших n становится неустойчивым.
Я попробовал в правой части использовать не $f_{n-1}$, а $f_n$ и получил
$f_n = f_{n-1}/(1 - a/n + \delta\cdot dt)$
Этот процесс оказался очень устойчивым, но ведет себя в высшей степени странно -
в зависимости от показателя $a$, он "развивается" то в положительной, то в отрицательной полуплоскости.
В последнем случае график оказывается просто отраженным относительно оси $t$.
Сама форма графика отлично передает поведение искомой функции.

Хотел бы услышать мнения о корректности такого метода.
Любопытно также то, что $f(0) = 0$, но процесс начать с нулевого значения не получится.
Я выбираю "разумные" значения $f_0$ типа $0.0001$, при этом меняется масштаб графика ( с сохранением пропорций).
Какое стартовое значение нужно брать, чтобы получать правильные значения функции?

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

[ewert]

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

$(f_n - f_{n-1})/dt = (a/t_n - \delta)\cdot f_{n-1}$
$f_n = (1 + a/n - \delta\cdot dt)\cdot f_{n-1}$
Такой процесс работает, но при больших n становится неустойчивым.

С какой стати ему быть неустойчивым, когда это всего-навсего метод Эйлера.

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

Я попробовал в правой части использовать не $f_{n-1}$, а $f_n$ и получил
$f_n = f_{n-1}/(1 - a/n + \delta\cdot dt)$
Этот процесс оказался очень устойчивым, но ведет себя в высшей степени странно

А это называется неявным методом Эйлера; он, может, чуток и устойчивее, но это непринципиально. Принципиально же то, что начальное условие нельзя задавать в нуле, т.к. это -- особая точка для дифференциального уравнения. Уж как минимум начинать в неявном методе можно не ранее $n>a$. В явном -- формально можно начать и с $n=1$, но фактически выйдет бред, поскольку на начальных шагах степенное слагаемое меняется на величины, много большие шага, так что о какой аппроксимации вообще можно говорить.

Но главное в другом -- всё это просто не нужно. Ибо:

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

Вычислять их прямо каждую $1/44100$ долю секунды слишком жирно.

Если и жирно, то с точностью до наоборот. Двадцати или пятидесяти тысяч тактов на вычисление и всего-то одного логарифма и одной экспоненты хватит за глаза, притом много раз. Это лет двадцать назад, может, и было актуально.

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

0. Полностью присоединяюсь к предыдущему оратору касательно излишества такого приёма для генерации звука на РС.
Время вычисления логарифма, как и экспоненты, примерно в 4 раз больше времени элементарной операции сложения или умножения (а для синуса/косинуса - приблизительно в 8). Экономия тут возможна на спецпроцессорах, где нет аппаратной реализации этих функций, или на очень старых PC, без сопроцессора. Однако в специализированных системах (например, электромузыкальный инструмент на 8-битном процессоре или аппаратура связи, на сигнальных процессорах) такое преобразование может быть актуально.
1. Общий подход - смотреть "исчисление конечных разностей". Где не делается допущение о пренебрежимой малости шага. Уравнения в конечных разностях имеют теорию, очень похожую на теорию дифуравнений. Оттуда и выводить схеиы вычислений.
2. Частный вопрос вычисления представленной функции.
Для PC совсем просто. Логарифм там достаточно быстр, как и экспонента.
А вот для процессора, у которого четыре действия, но деление слишком медленное... (реально для некоторых сигнальников, у которых все операции в 1 такт а деление в 16).
Я бы отдельно считал первый, отдельно второй сомножители. Второй - просто умножением на каждом шаге на константу.
А для первого пришлось бы приближённо логарифмировать. Примерно как в a-law или mu-law в телефонии. Сдвигом - двоичный показатель и линейным приближением для мантиссы. И потом приближённо потенцировать.

[Roman Voznyuk]

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

Если и жирно, то с точностью до наоборот. Двадцати или пятидесяти тысяч тактов на вычисление и всего-то одного логарифма и одной экспоненты хватит за глаза, притом много раз. Это лет двадцать назад, может, и было актуально.

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

[profrotter]

Lesobrod у вас параметр $a$ натуральное число?

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

Ну, для натурального $a$ совсем уж тривиально...