Еще раз спрашиваю Горина и Ost - вы оба уверены, что решаете именно задачу, где источник движется в среде распространения волн, а приемник неподвижен ?
Если уверены в этом, то смею вас обоих огорчить, что ваше "точное" решение - НЕ точное для этого случая.
Неточность началась отсюда:
\[ L{i+1}=\sqrt{L_i^2+v^2T_0^2-2L_ivT_0 fCos_i} \]
\[ fCos_{i+1}=\frac{L_i fCos_i - vT_0}{L_{i+1}} \]
Это вы рассчитываете ТЕКУЩЕЕ геометрическое расстояние между источником и приемником/
Вернее, если перейти к разности \[ L_{i+1}-L_i \], то это изменение текущего расстояния между источником и приемником за время одного собственного периода источника
А для решения ВАШЕЙ задачи нужно определять ЗАПАЗДЫВАЮЩЕЕ расстояние до источника на момент приема, и соответственно разность ЗАПАЗДЫВАЮЩИХ расстояний.
Даже если вы возразите, что имели в виду Lo как исходное запаздывающее расстояние, то все равно формула останется НЕ точной, потому что запаздывающее расстояние в ВАШЕЙ задаче будет изменяться НЕ со скоростью v (на вектор перемещения
vTo), а ПО ДРУГОМУ (разница с вашим решением второго порядка по скорости).
Для дискретной (пошаговой) функции угла замечание аналогичное:
проекция запаздывающего расстояния L на ось Х (по условию параллельной движению источника) НЕ изменяется в точности на перемещение vTо - изменение ДРУГОЕ.
Проще всего эта разница видна для случая, когда источник стартует от приемника в момент to=0 со скоростью v.
Приемник в последующие текущие моменты времени t будет воспринимать сигналы от источника с запаздывающего расстояния L', для которого можно записать соотношение
\[ L'=vt-vL'/c \]
откуда
\[ L'=\frac{vt}{1+v/c} \approx vt(1-v/c) \]
Явно видно, что за какое-то время dt запаздывающее расстояние изменяется НЕ на величину перемещения vdt, а сложнее.
Если приемник НЕ на линии движения источника, то зависимость еще сложнее и тоже НЕ сводится к тому что вы оба получили.
Ваши же сравнительно простые запаздывающие решения являются точными для другого случая - когда источник неподвижен в среде распространия волн, а движется приемник - тогда действительно запаздывающее расстояние для приемника будет изменяться на величину его перемещения vdt.