Система связанных уравнений в частных производных в атомной физике

Aug 19 2020

Мой вопрос касается реализации системы связанных PDE с подпрограммой Mathematicas NDSolve. Я рассматриваю одномерную игрушечную модель в атомной физике. Модель описывает два поля$\psi =\psi(t,z)$ и $\sigma= \sigma(z;t)$ соединены друг с другом, т.е. $$ i \hbar \partial_t \psi = -\frac{\hbar^2 }{2 m} \psi_{zz} +V \psi +\frac{\hbar^2 \alpha_s }{m}\sigma^{-2} \left| \psi \right|^2 \psi+\frac{\hbar^2}{2m }\sigma^{-2}\psi+\frac{1}{2} m \omega_{\perp} \sigma^2 \psi +\frac{\hbar^2 }{2 m} \sigma^{-2}\sigma_z^2 \psi , \\ 0 =-\frac{\hbar^2}{4 m}\sigma \sigma_{zz}+\frac{\hbar^2 }{ m } \sigma^{-3} \sigma_z^2 -\frac{\hbar^2 }{4 m} \sigma \sigma_z \frac{1}{\left| \psi \right|^2} \left(\psi\psi_z^*+\psi^* \psi_z\right)+\frac{\hbar^2}{2 m }\sigma^{-3}-\frac{m \omega_{\perp}}{2} \sigma + 2 \frac{\hbar^2 \alpha_s}{m } \sigma^{-3} \left| \psi \right|^2 $$ Дополнительно я накладываю периодические граничные условия для $\psi(-L/2,t) = \psi(L/2,t)$ и $\sigma(-L/2,t) = \sigma(L/2,t)$ и зададим некоторые начальные условия $\psi(z,0)=f(z)$ и $\sigma(z,0)=g(z)$.

ИЗМЕНЕНО:

Вот моя текущая версия кода

    (*constants*)
h = 1; (* Planck constant *)
m = 1; (* particle mass *)
Subscript[\[Alpha], s] = 1; (* scattering length *)
\[Omega] = 1; (* frequency *)
V = 0; (* potential *)

(*ranges*)
L = 2; (*length of the box *)
tmin = 0;
tmax = 0.1;

(*equations*)
eqn1 = I  D[\[Psi][z, t], t] == -h^2/(2 m) D[\[Psi][z, t], z, z] + 
    V \[Psi][z, t] + 
    h^2 Subscript[\[Alpha], s]/
      m  \[Sigma][z, t]^(-2) Abs[\[Psi][z, t]]^2 \[Psi][z, t] + 
    h^2/(2 m) \[Sigma][z, t]^(-2) \[Psi][z, t] + 
    m \[Omega] /2 \[Sigma][z, t]^2 \[Psi][z, t] + 
    h^2/(2 m) \[Sigma][z, t]^(-2) D[\[Sigma][z, t], z]^2 \[Psi][z, t];

eqn2 = -h^2/(4 m) \[Sigma][z, t]  D[\[Sigma][z, t], z, z] ==  
   h^2/(2 m) \[Sigma][z, t]^(-3) D[\[Sigma][z, t], z]^2 -  
    h^2/(4 m) \[Sigma][z, t]   D[\[Sigma][z, t], z]  /
      Abs[\[Psi][z, t]]^2  ( \[Psi][z, t]  D[\[Psi][z, t], 
         z] + \[Psi][z, t] D[\[Psi][z, t], z]) + 
    h^2/(2 m) \[Sigma][z, t]^(-3)   - m \[Omega] /2 \[Sigma][z, t] + 
    2 h^2 Subscript[\[Alpha], s]/
      m \[Sigma][z, t]^(-3) Abs[\[Psi][z, t]]^2;

(*boundary conditions*)
bc = \[Psi][L/2, t] == \[Psi][-L/2, t];
bcwidth = \[Sigma][L/2, t] == \[Sigma][-L/2, t];

(*initial conditions*)
icwidth = \[Sigma][z, 0] == z^2 + 1;
icdwidth = D[\[Sigma][z, t], t] == 2 /. t -> 0;
icwave = \[Psi][z, 0] == Exp[-((z)^2)];

(*solve system*)
sol1 = NDSolve[{eqn1, eqn2, bc, bcwidth , icwave, icwidth, 
    icdwidth}, {\[Psi], \[Sigma]}, {z, -L/2, L/2}, {t, tmin, tmax}, 
   Compiled -> True, MaxSteps -> {500, Infinity}];

К сожалению, это связано с двумя проблемами, первая касается самого решателя, поскольку в моем уравнении для второго поля нет производной по времени. $\sigma$ он обрабатывает систему как DAE и выдает два предупреждения

NDSolve :: pdord: Некоторые функции имеют нулевой дифференциальный порядок, поэтому уравнения будут решаться как система дифференциально-алгебраических уравнений. >>

NDSolve :: mconly: Для метода IDA доступен только машинный реальный код. Невозможно продолжить со сложными значениями или исключениями с плавающей запятой. >>

Я не знаю, настоящая ли это проблема (я использую Mathematica 9.x). Второй более проблематичный, касается количества используемых точек сетки. Я думаю, это в основном исходит из самих уравнений и вызывает ошибку, заключающуюся в том, что он не может найти подходящее решение в пределах допуска.

NDSolve :: mxsst: Использование максимального количества 500 точек сетки, разрешенных параметрами MaxPoints или MinStepSize для независимой переменной z. >>

NDSolve :: icfail: Не удалось найти начальные условия, которые удовлетворяют функции невязки в пределах указанных допусков. Попробуйте задать начальные условия как для значений, так и для производных функций. >>

Я также попытался дать ему дополнительные исходные данные, как указано в сообщении об ошибке, но безуспешно. Вопрос. Я не знаю, есть ли какой-либо потенциал для улучшения моего кода, или обновление до новой версии Mathematica решит проблему, или, в худшем случае, это «слишком уродливая» система для числовой обработки.

Ответы

1 AlexTrounev Aug 26 2020 at 21:54

Для решения такого рода задач волновую функцию можно разделить на две части. $\psi=\psi_1+i\psi_2$. Также мы используем некоторые варианты, NDSolveпозволяющие решить эту проблему. Предположим, что$\sigma$ реально, тогда у нас есть