原子物理学における結合偏微分方程式システム

Aug 19 2020

私の質問は、MathematicasNDSolveルーチンに結合されたPDEのシステムの実装についてです。原子物理学の1次元トイモデルを考えています。モデルは2つのフィールドを記述します$\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}];

残念ながら、2つの問題があります。最初の問題は、2番目のフィールドの方程式に時間微分がないため、ソルバー自体に関係します。 $\sigma$ システムをDAEとして処理し、この2つの警告を出します

NDSolve :: pdord：一部の関数は微分次数がゼロであるため、方程式は微分代数方程式のシステムとして解かれます。>>

NDSolve :: mconly：メソッドIDAの場合、マシンの実際のコードのみが使用可能です。複素数値または浮動小数点例外を超えて続行できません。>>

これが「本当の」問題かどうかはわかりません（私はMathematica9.xを使用しています）。2つ目はもっと問題があり、使用されるグリッドポイントの量に関係します。これは主に私が推測する方程式自体に由来し、許容範囲内で適切な解を見つけることができないというエラーを引き起こします。

NDSolve :: mxsst：独立変数zのMaxPointsまたはMinStepSizeオプションで許可されているグリッドポイント500の最大数を使用します。>>

NDSolve :: icfail：指定された許容範囲内で残差関数を満たす初期条件を見つけることができません。関数の値と導関数の両方の初期条件を与えてみてください。>>

また、エラーメッセージで示されているように、追加の初期データを提供しようとしましたが、成功しませんでした。質問私が知らないのは、コードを改善する可能性があるのか、それとも新しいバージョンのMathematicaにアップグレードすることで問題が解決するのか、最悪の場合、数値処理には「醜い」システムなのかということです。

回答

1 AlexTrounev Aug 26 2020 at 21:54

この種の問題を解決するために、波動関数を2つの部分に分割することができます $\psi=\psi_1+i\psi_2$。またNDSolve、この問題を解決可能にするためにいくつかのオプションを使用します。それを仮定しましょう$\sigma$ 本物であるなら、