Радиационное трение

[AlexFuser]

Здравствуйте! Существует известное уравнение движения частицы, учитывающее радиационное трение.
$m\ddot{r(t)}=F_{rad}+F_{ext}(t)=\frac{q^2}{6\pi\varepsilon_{0}c^3}\dddot{r(t)}+F_{ext}(t)$
И если считать, что $F_{ext}=0$, то возникают так называемые "runaway solutions"
$\ddot{r}=\tau\dddot{r}\ \Rightarrow r(t)=r_{0}+v_{0}t+a_{0}e^\frac{t}{\tau}$
А если, допустим, поместить частицу в постоянное электромагнитное поле, т. е $F_{ext}=q[E+v(t)\times B]$,
$m\ddot{r(t)}=\frac{q^2}{6\pi\varepsilon_{0}c^3}\dddot{r(t)}+q[E+v(t)\times B]$
то кажется, что движение частицы по спирали будет затухать.
Я хотел получить графики зависимостей: $v(t), a(t)$
Однако, численно решая соответствующие диффуры(опыта у меня в этом немного), я все равно получаю бесконечные runaway решения, а не ожидаемое затухание. Возможно, моя ошибка заключается в неверном подборе начальных условий. Прошу помощи, если понимаете какие начальные условия можно принять. Возможно,у меня есть ошибки в рассуждениях

[sergey zhukov]

AlexFuser
По моему, "runaway solution" - это просто обращенное во времени, нефизичное по сути, решение, т.к. сила радиационного трения тут взята с противоположным знаком (заряд не излучает и тормозится, а поглощает излучение и ускоряется). Такое же решение, как и обращение во времени решения для механической системы с трением. Следует рассматривать уравнение вида: $\ddot r=-\tau\dddot r$
Аналогично и для второго случая.

[mihiv]

Сила радиационного трения мала по сравнению с силой Лоренца, поэтому, мне кажется , применим метод последовательных приближений. В нулевом приближении считаем $F_{rad}=0$ и находим траекторию движения заряда $\vec r_0(t)$. Затем в уравнении 1-ого приближения в выражение для $F_{rad}$ подставляем $ \dddot r_0(t) $ и т.д.

[revos]

Непростое это дело.
Академик пишет https://www.mathnet.ru/links/3f2bd92151 ... n11203.pdf

[zykov]

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

[svv]

Попробуйте, исходя из дифференциального уравнения и конечных условий $\mathbf r(t_0),\dot{\mathbf r}(t_0), \ddot{\mathbf r}(t_0)$, строить решение от момента $t_0$ назад во времени. Понятно, что это несколько другая задача, и на какие-то Ваши вопросы она не ответит; но на какие-то всё же ответит.
Вопрос, откуда взять конечное ускорение $\ddot{\mathbf r}(t_0)$, не так важен — хотя бы из условия $\dddot{\mathbf r}(t_0)=0$. Ошибка от неправильного задания ускорения будет сказываться только в окрестности $t_0$, в том смысле, что если эту окрестность отбросить, остальная часть решения будет одним из правильных движений частицы в присутствии силы радиационного трения.

[LLeonid3]

AlexFuser насколько помнится, радиус кривизны (до релятвиской поправки) не зависит от скорости, но только от магнитного поля, заряда и массы частицы, т. е. заряженная частица будет двигаться по одной окружности, а с потерей энергии с разным периодом. Поэтому считать надо по углу поворота, а время просто суммировать, пока частица не достигнет минимальной скорости, предположим $1.0$ мм/сек

[svv]

Кстати, случай, когда поле $\mathbf E=0$, поле $\mathbf B$ однородно и частица движется в плоскости, перпендикулярной $\mathbf B$, решается аналитически.

[sergey zhukov]

LLeonid3
Сила Лоренца пропорциональна скорости, а центростремительное ускорение при фиксированном радиусе - квадрату скорости. Так что с ростом скорости радиус окружности возрастает.

[LLeonid3]

sergey zhukov да у меня заскок случился, думал о периоде, а написал радиус
Это период вращения заряженной частицы НЕ зависит от её скорости, а радиус "орбиты" как раз и меняется от скорости частицы (энергии), поэтому и числено интегрировать надо по углу, а не по времени, погрешность получится на порядки меньше.
Программку что-ли сделать

[Cos(x-pi/2)]

Продолжение из этой темы. Рассмотрим частный пример движения частицы с зарядом $q$ и массой $m:$ пусть $\mathbf{E}=0,$ частица движется в плоскости, перпендикулярной вектору постоянного однородного магнитного поля $\mathbf{B}.$

Орт декартовой системы координат $\mathbf{e}_z$ выберем вдоль $\mathbf{B}.$ Тогда $\mathbf{B}=B\mathbf{e}_z,$ радиус-вектор частицы есть $\mathbf{r}(t) = x(t)\mathbf{e}_x+y(t)\mathbf{e}_y.$

Уравнение Ньютона (в гауссовой системе единиц)

$$m \ddot{\mathbf{r}} = \dfrac{q}{c}[\dot{\mathbf{r}} \times \mathbf{B}]+\dfrac{2q^2}{3c^3} \dddot{\mathbf{r}}$$

сводится к системе двух уравнений для координат $x(t)$ и $y(t)$ частицы:

$$\ddot{x}=\omega_c \dot{y} + \tau \dddot{x}$$

$$\ddot{y}=-\omega_c \dot{x} + \tau \dddot{y},$$

где

$$\omega_c = \dfrac{qB}{mc}\,, \qquad \tau = \dfrac{2q^2}{3mc^3}\,.$$

При $\tau=0$ частица вращалась бы по окружности постоянного радиуса с "циклотронной частотой" $\omega_c.$ Четыре постоянных интегрирования (начальные данные): координаты центра орбиты $x_c$ и $y_c$, начальная фаза вращения $\alpha$, радиус орбиты $C.$ Полагая для удобства, что начало координат помещено в центр орбиты (так что $x_c=y_c=0),$ можно один раз проинтегрировать уравнения с равными нулю постоянными интегрирования; такие же уравнения получатся для переменных $x-x_c$ и $y-y_c$ с произвольными $x_c,\,y_c.$

Поступим так же при $\tau > 0.$ Вот получившиеся дифференциальные уравнения (ДУ) в записи с безразмерной переменной времени $t'=\omega_c t$ и с безразмерным неотрицательным параметром радиационного трения $w=\omega_c \tau$

$$\dot{x} = y + w\ddot{x},$$

$$\dot{y} = -x + w\ddot{y}.$$

1. Сначала найдём точное решение, чтобы было с чем сравнивать приближённые ответы. "Методом комплексных амплитуд" ищем решение в виде

$$x = \operatorname{Re}\left (A\,\exp(i\Omega t')\right ),$$

$$y = \operatorname{Re}\left (C\,\exp(i\Omega t')\right ),$$

где комплексные постоянные $A,\,C,\,\Omega$ определяются системой уравнений

$$(\Omega^2w+i\Omega)A-C=0$$

$$A+(\Omega^2w+i\Omega)C=0.$$

Такая система имеет не равные нулю решения $A,C$ только тогда, когда её определитель равен нулю:

$$(\Omega^2w+i\Omega)^2+1=0$$

Решение этого уравнения для $\Omega$ можно выбрать так, что:

$$\Omega^2w+i\Omega=i,\qquad A=-iC,$$

$$\Omega = \dfrac{i}{2w}\left ( -1 +\sqrt{1-i4w} \right )$$

Обозначим: $G$ - мнимая часть, $R$ - вещественная часть частоты $\Omega,$ т.е.

$$\Omega = R + iG$$

Положим $C=1$ - это означает, что начальную фазу $\alpha$ выбираем равной нулю, а расстояние от начала координат до частицы при $t=0$ выбираем за единицу длины. Тогда $A=-i.$

Таким образом, вот точные выражения для траектории частицы c $\alpha=0,\,C=1:$

$$x(t')=e^{-Gt'}\sin(Rt'),$$

$$y(t')=e^{-Gt'}\cos(Rt').$$

Для перехода к исходным единицам достаточно умножить $x$ и $y$ на произвольный вещественный коэффициет $C$ с размерностью длины и подставить $t'=\omega_c t,$ где $\omega_c=\frac{qB}{mc}.$ К аргументам синуса и косинуса можно добавить произвольное вещественное число $\alpha.$ Если к получившимся $x(t)$ и $y(t)$ добавить произвольные значения координат центра орбиты $x_c$ и $y_c,$ то в ДУ следует под координатами $x$ и $y$ понимать разности $x-x_c$ и $y-y_c.$

В решении $\Omega = R+iG$ уравнения $\Omega^2w+i\Omega=i$ величина $R$ имеет смысл частоты вращения, а $G$ - параметр "затухания спирали"; оба параметра безразмерные, измеряются в долях от $\omega_c.$ Из уравнения для $\Omega$ следует, что $w(R^2-G^2)-G=0$ и $2wRG+R=1.$ Видно, что в отсутствие трения, т.е. при $w=0$ получается, как и должно быть, $R=1,\,G=0,$ а при небольшом трении, т.е. при $0<w\ll 1$ получается немного меньшее единицы значение $R$ и маленькое положительное $G.$ Значит, с трением спираль будет не раскручиваться, а сжиматься, "затухать".



2. Теперь для упражнения попытаемся решать ту же систему ДУ приближённо.

Оказывается, простой метод конечных разностей ("метод Эйлера") даже в отсутствие трения плохо работает: вместо стационарной окружности у меня получалась постепенно раскручивающаяся спираль. Поэтому применяю более совершенный способ приближённых вычислений, "метод Рунге-Кутта". А именно:

Переходим к системе ДУ первого порядка c вспомогательными функциями $p(t')$ и $r(t')$:

$$\dot{x}=p,\quad \dot{p}=(p-y)/w, \quad \dot{y}=r, \quad \dot{r}=(r+x)/w.$$

Шаг по безразмерному времени $t'$ обозначаю как $h.$ Приращения функций записываю в виде оборванного ряда по степеням $h$ с точностью до членов порядка $h^4$ включительно; например, для функции $p$ имеем приближённое равенство:

$$p(t'+h)=p(t')+\dot{p}h+\frac{1}{2}\ddot{p}h^2+\frac{1}{6}\dddot{p}h^3+\frac{1}{24}\ddddot{p}h^4,$$

где производные относятся к точке $t'.$ Затем эти производные последовательно выражаем через сами функции в точке $t'$ с помощью наших ДУ; например, $\dot{p}=(p-y)/w,$ дифференцируя это ещё раз и снова пользуясь ДУ, получаем выражение для $\ddot p,$ и так далее.

Значения функций в точке $t'$ пишем с индексом, например, $n.$ В точке $t'+h$ они тогда будут с индексом $n+1.$ То есть: $p_{n+1}=p_n +...,$ и так далее. Так получается система рекуррентных соотношений, позволяющих в цикле по $n$ от нуля до какого-нибудь большого $N\gg 1$ вычислять интересующее нас множество значений координат $x_n,\,y_n,$ вместе со вспомогательными $p_n, \, r_n,$ имеющими смысл значений скорости $\dot{x}_n,\,\dot{y}_n.$

a) Посмотрим для начала, что получается при положительном шаге $h=0.01$ с параметром трения $w=0.2.$

Пусть начальные координаты частицы: $x_0=0,\, y_0=1.$ При этом из приведённых выше формул точного решения (c $\alpha=0,\, C=1)$ видно, что начальную скорость можно выбрать так: $p_0=\dot x(0) = R,$ $r_0=\dot y(0) = -G,$ где $R$ и $G$ определяются из решения уравнения для $\Omega=R+iG.$ Вычислить $\Omega$ поручаю Маткаду (но не знаю, насколько хорошо он с таким заданием справляется :). Вот маткадные фрагменты расчёта, рекуррентные формулы:



Ниже на графиках красная линия - точное решение, пунктир - результат приближённого вычисления по указанным рекуррентным формулам. Пунктир это способ изображения графика, а не сами вычисленные точки. Приближённых точек так же много, как и точных ($N=2000),$ но если их все изобразить, то они полностью закрашивают собой предыдущий график. Результат при $h=0.01$ и $w=0.2:$



Видно, что начало приближённого решения соответствует точному решению (пунктир идёт по красной линии), но затем быстро развивается неустойчивость приближённого расчёта. По физическому смыслу величина $w$ должна быть малой. Но при меньших значениях $w$ картина становится только хуже: неустойчивость проявляется ещё раньше.

б) Выполняем расчёт "назад во времени" - при отрицательном $h=-0.01,$ с прежним значением параметра трения $w=0.2.$

Начальные значения координат и скорости в этом примере выбираю те, которые получаются из точного решения для конечного момента времени $t'=20,$ т.е. для точки с номером $N=2000:$ $x_0=-4.0379830047\cdot 10^{-3},$ $y_0=3.3374190164\cdot 10^{-2},$ $p_0=3.1938495138\cdot 10^{-2},$ $r_0=-1.88\cdot 10^{-3}.$ После завершения цикла вычислений (с прежними рекуррентными формулами, по $n$ от нуля до $N=2000,$ но теперь с отрицательным шагом $h<0),$ все получившиеся координаты $x_n,\, y_n$ перенумеровываем в обратном порядке. И... ура, неустойчивость исчезла, результат соответствует точному решению:



"Затухающая спираль" на плоскости $x,y$ в этом примере выглядит так:



в) Расчёт "назад во времени", при $h=-0.01,$ с прежним значением параметра трения $w=0.2,$ но с произвольными начальными условиями:

Начальные условия теперь задаю как попало; например, пусть $x_0=y_0=p_0=r_0=0.1.$ По отношению к нумерации "вперёд во времени" эти условия конечные - они для точки $t'$ с номером $N=2000.$ После завершения такого же, как и раньше, цикла вычислений по $n$ от нуля до $N=2000,$ с отрицательным шагом $h=-0.01,$ получившиеся координаты $x_n,\, y_n$ перенумеровываем в обратном порядке. В результате, получившиеся $x_0,\,y_0$ относятся уже к точке $t'=0.$ По ним определяем начальный радиус $C=\sqrt{x_0^2+y_0^2}$ и начальную фазу $\alpha=\arctg \frac{x_0}{y_0}$ для сравнения с точным решением

$$x(t')=Ce^{-Gt'}\sin(Rt'+\alpha),$$

$$y(t')=Ce^{-Gt'}\cos(Rt'+\alpha).$$

Оказывается, в этом случае результат приближённого вычисления тоже соответствует точному решению:



Лишь на первых шагах расчёта "назад во времени" (т.е. вблизи $t'=20$ в этом примере) заметно обусловленное произвольными начальными условиями большое отличие приближённого результата от точного решения:




г) Расчёт "назад во времени" при $h=-0.01$ с малым значением параметра трения: $w=0.01.$

Пусть максимальное $t'$ будет равно $1000,$ так что $N=1000/|h|=100000.$ Начальные условия (т.е. условия для точки $t'$ с номером $N)$ задаю произвольно: $x_0=y_0=p_0=r_0=0.00033.$ Сравнение с точным решением делается тем же способом, что и в предыдущем примере; только теперь $N$ в пятьдесят раз больше. Цель этого теста - убедиться, что даже на длинном куске траектории приближённый расчёт не сбивается, согласуется с точным решением. Оказывается, так и есть. Вот часть траектории в начале временного интервала:



Последняя часть траектории:



Вся "спираль" на плоскости $x,y$ в этом примере выглядит так:



Маткад сумел даже кое-как анимировать примерно 80 оборотов по этой траектории (здесь ссылка на avi-файл, некоторые браузеры могут его показывать, но не знаю, все ли ).


О перенумеровке точек траектории в обратном порядке:

Есть варианты. Вычисленные по указанным в начале рекуррентным формулам (c $h<0)$ массивы $x_n, y_n$ можно в цикле по $j$ от нуля до $N$ переписать в новые массивы, например, $x1_j:=x_{N-j}, \, y1_j:=y_{N-j},$ а затем, чтобы вернуться к исходным обозначениям, в отдельном цикле записываем их на место исходных величин: $x_j:=x1_j$ и т.д. Более простой вариант: в рекуррентных формулах можно дать массивам новые обозначения, например, $x1_n,\,y1_n$ вместо $x_n, y_n$ и т.д.; тогда для перенумеровки в обратном порядке достаточно в цикле по $j$ от нуля до $N$ записать $x_j:=x1_{N-j},\,y_j:=y1_{N-j}.$ Наконец, самый простой вариант - вести рекуррентные вычисления сразу в обратном порядке, т.е. в цикле по $n$ от $N$ до $1$ с шагом $-1$ вычислять $x_{n-1}=\text{выражение с индексом n}$ и т.д.; тогда не надо будет ничего перенумеровывать.

[svv]

Cos(x-pi/2)
Спасибо за Ваш труд.

[Cos(x-pi/2)]

svv
Огромное спасибо Вам!

Мне не за что (я лишь написал про то, что попытался понять на своём ученическом уровне из Ваших сообщений ).

Все идеи - в Вашем потрясающем сообщении:

svv в сообщении #1588186 писал(а):

Попробуйте, исходя из дифференциального уравнения и конечных условий $\mathbf r(t_0),\dot{\mathbf r}(t_0), \ddot{\mathbf r}(t_0)$, строить решение от момента $t_0$ назад во времени. Понятно, что это несколько другая задача, и на какие-то Ваши вопросы она не ответит; но на какие-то всё же ответит.

Вопрос, откуда взять конечное ускорение $\ddot{\mathbf r}(t_0)$, не так важен — хотя бы из условия $\dddot{\mathbf r}(t_0)=0$. Ошибка от неправильного задания ускорения будет сказываться только в окрестности $t_0$, в том смысле, что если эту окрестность отбросить, остальная часть решения будет одним из правильных движений частицы в присутствии силы радиационного трения.

[AlexFuser]

Cos(x-pi/2), я разобрался, спасибо! Я Вам очень благодарен!
svv, огромное спасибо!
Для меня это было очень полезно! Еще раз большое спасибо.

[svv]

AlexFuser, Вам тоже спасибо за интересную задачу!

(Оффтоп)

Думал, не изложить ли своё точное решение с ложечкой дёгтя, но не уверен, что это нужно.