дифференциальное уравнение с большими коэффициентами

[Anatoly89]

Всем здравствуйте!

Пытаюсь решить в матлабе дифур вида $\ddot{\psi}+\frac{1}{\tau}\dot{\psi}+\Omega^2\sin \psi=0,~1/\tau=1.2108\times10^7,~\Omega=2.5217\times 10^8$

Использовал ф-ции ode45 и ode23s. Первая функция зависает при таких болиших коэффициентах, а вторая выдает ошибку

In linThiele_solve at 3
Warning: Matrix is close to singular or
badly scaled.
Results may be inaccurate.
RCOND = 1.000000e-16.

Что делать в таких случаях?

[Klad33]

Эх, научиться бы решать в общем виде:

y''+a y' +b sin(y)=0

Ни один пакет мне не дал ничего путного. Даже Вольфрам_Альфа
Сделал запрос на форум exponenta.ru . Может, там дадут дельный совет. Но не раньше завтра

[AlexValk]

Начальные/граничные условия у Вас есть или нужны свойства общего решения?
На каком интервале (в какой области $t$) Вы ищите решение?
Есть ли представление о том, в каких пределах лежит $\psi(t)$?

Математическое обсуждение:
Для ДУ с большими/малыми параметрами существует своя наука (асимптотические методы), позволяющая получать приближенные аналитические решения. Вообще-то, ДУ с большими/малыми параметрами обычно трудная задача для любых универсальных программ численного решения ДУ. Проблему частично можно почувствовать, рассмотрев такое, тривиальное с аналитической точки зрения уравнение, как $y'(t)=\sin(10^{12} t), y(0)=1, t\in[0;1]$ - сколько периодов синуса должна "прокопать носом" программа, прежде чем она добереться от точки 0 до точки 1?!
Однако, используя идеи асимптотической теории численным программам часто можно помочь, предварительно преобразовав ДУ к виду, более пригодному для численных расчетов. К сожалению, подбор таких аналитических преобразований часто неоднозначен и неформален - здесь требуется опыт, а иногда и искусство.
Например, в Вашем случае $\Omega$ и $\tau^{-1}$ - большие величины одного порядка, и получается, что коэффициенты при производных в ДУ с ростом порядка убывают примерно пропорционально большому параметру. Поэтому здесь будет эффективна процедура ``сжатия времени'' - замена переменной $t=\tau T$. Для функции $\Psi(T)=\psi(t)$ получится ДУ

$$\Psi''+\Psi'+\alpha^2\sin(\Psi)=0,$$

где $\alpha=\Omega \tau$ - параметр, который в Вашем случае ($\alpha\approx 20$) можно считать "нормальным'' (не очень большим и не маленьким). Тогда:
1. Для вычисления решения в области $t\sim \tau$, можно использовать программы типа упоминаемых Вами ode45 - решение ДУ с таким параметром $\alpha$ не должно вызывать у численных алгоритмов трудностей, вплоть до значений $t$ для которых $t/\tau$ будет несколько десятков и более. Проблема здесь для хорошей программы может состоять только в том, что ей (программе) нужно будет сделать слишком много шагов, чтобы дойти до нужной нам очень далеко расположенной точки $t$ (и при этом убедиться, что в этой точке решение практически равно 0 - см след. пункт.).
2. Для вычисления (фактически - построения) решения в области $t\gg \tau$ можно использовать аналитические асимптотические методы, основанные на замене нелинейного уравнения $\Psi''+\Psi'+\alpha^2\sin(\Psi)=0$ линейным $\Psi''+\Psi'+\alpha^2(\Psi-2\pi n)=0$. В частности, для $\alpha>1/2$ главный член асимптотики решения при больших $t/\tau$ имеет экспоненциально-осциллирующий характер

$$\psi(t)\sim 2\pi n+A_0 e^{-\frac{t}{2\tau}}\cos\left(\sqrt{\Omega^2-\frac{1}{4\tau^2}}\, t+\varphi_0\right).$$

Целое число ("индекс") $n$, "амплитуда" $A_0$ и "фазовый сдвиг" $\varphi_0$ определяются начальными условиями. Эти параметры можно найти "сшивая'' численные решения, полученные в области $t\sim \tau$ с аналитической асимптотической формулой из области $t\gg \tau$. В качестве точки сшивки можно использовать $t=15\tau$ если достаточна точность $\sim 10^{-4}$ или большие значения $t\approx -2\ln(\varepsilon)\tau $ - если нужна более высокая точность $\varepsilon$.
При $\alpha\leq 1/2$ асимптотический вид решения имеет чисто экспоненциально-убывающий характер:

$$\psi(t)\sim 2\pi n+A_0 e^{-\frac{t}{2\tau}}\exp\left(-\sqrt{\frac{1}{4\tau^2}-\Omega^2}\, t\right).$$

Эта ситуация ($t\gg \tau$) может представлять интерес в области только при $\alpha\approx 1/2$, иначе во второй экспоненте стоит огромный отрицательный показатель и решение практически неотличимо от нуля.
3. Интересным моментом в задаче являются бифуркации - скачкообразные изменения целого числа $n$ на 1 при непрерывном изменении начальных условий.

Физическое обсуждение:
По внешнему виду данная задача - это классическое уравнение математического маятника (без предположения о малости амплитуды колебаний угла отклонения $\psi$) с добавочным учетом простейшего варианта трения (сила трения предполагается пропорциональной скорости). Поэтому многие свойства решения можно понять (и говорить о них) на физическом языке.
A. Начальное условие $\psi(0)=\psi_0$ можно рассматривать, без ограничения общности, только в интервале $\psi\in[0;2\pi]$ - в противном случае достаточно сделать замену $\psi(t)\to \psi(t)-2\pi n$, где $n=[\psi_0/2\pi]$ - целое число.
B. Величина $\Omega/2\pi$ имеет смысл частоты колебаний (при малости амплитуды и в отсутствие трения), а $\tau$ - характерного времени затухания.
C. Решение может быть интересно только на временах, когда $t$ не превышает характерное время затухания $\tau$ слишком значительно (иначе $\psi$ практически равно 0).
D. Слагаемое $2\pi n$ в асимптотике решения - при значительной начальной скорости маятника (а также если скорость не так велика, но исходное положение маятника далеко от равновесия) он совершит несколько (целое число $n$) оборотов, прежде чем “успокоится”.
E. Осцилляции асимптотики решения при $\alpha>1/2$ – если трение не слишком велико, то после нескольких оборотов маятник за счет трения "успокоиться" и начинает колебаться вокруг точки равновесия, постепенно уменьшая амплитуду колебаний со временем - стремясь в итоге к положению равновесия.
F. Затухающий характер асимптотики решения при $\alpha>1/2$ – в случае достаточно сильного трения маятник будет завершать последний оборот в режиме "торможения" - приближаясь к положению равновесия с одной стороны (причем возможны оба случая: как слева, так и справа).
G. Бифуркации решения при плавном изменении начальных условий - плавно увеличивая начальную скорость можно добиться, чтобы маятник с трением сделал произвольное число оборотов, прежде чем он "успокоится". Это число оборотов может также на 1 измениться скачком при плавном изменении начального положения маятника, если его начальная скорость фиксирована.

[AlexValk]

Уточнение.
Чисто экспоненциальное поведение асимптотики решения относится к случаю $\alpha<1/2$.
А при $\alpha=1/2$ асимптотический вид решения имеет особый экспоненциально-степенной характер:

$$\psi(t)\sim 2\pi n+(A_0+B_0t) e^{-\frac{t}{2\tau}}.$$

(Кажется этот случай для линейного осциллятора с затуханием называется критическим).