Gekoppeltes PDE-System in der Atomphysik

Aug 19 2020

Meine Frage betrifft die Implementierung eines Systems gekoppelter PDEs an die NDSolve-Routine von Mathematicas. Ich betrachte ein eindimensionales Spielzeugmodell in der Atomphysik. Das Modell beschreibt zwei Felder$\psi =\psi(t,z)$ und $\sigma= \sigma(z;t)$ miteinander gekoppelt dh $$ 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 $$ Zusätzlich lege ich periodische Randbedingungen für fest $\psi(-L/2,t) = \psi(L/2,t)$ und $\sigma(-L/2,t) = \sigma(L/2,t)$ und stellen Sie einige Anfangsbedingungen ein $\psi(z,0)=f(z)$ und $\sigma(z,0)=g(z)$.

EDITIERT:

Hier ist meine aktuelle Version des Codes

    (*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}];

Leider gibt es zwei Probleme, das erste betrifft den Solver selbst, da meine Gleichung für das zweite Feld keine Zeitableitung enthält $\sigma$ Es behandelt das System als DAE und gibt diese beiden Warnungen aus

NDSolve :: pdord: Einige der Funktionen haben eine Differentialordnung von Null, daher werden die Gleichungen als System differenziell-algebraischer Gleichungen gelöst. >>

NDSolve :: mconly: Für die Methode IDA ist nur realer Maschinencode verfügbar. Komplexe Werte oder Gleitkomma-Ausnahmen können nicht fortgesetzt werden. >>

Ich weiß nicht, ob dies ein "echtes" Problem ist (ich verwende Mathematica 9.x). Der zweite ist problematischer, er betrifft die Anzahl der verwendeten Gitterpunkte. Dies kommt hauptsächlich von den Gleichungen selbst, die einen Fehler verursachen, dass er innerhalb der Toleranzgrenzen keine geeignete Lösung finden kann.

NDSolve :: mxsst: Verwenden der maximalen Anzahl von Rasterpunkten 500, die von den Optionen MaxPoints oder MinStepSize für die unabhängige Variable z zugelassen werden. >>

NDSolve :: icfail: Anfangsbedingungen, die die Restfunktion innerhalb der angegebenen Toleranzen erfüllen, können nicht gefunden werden. Versuchen Sie, Anfangsbedingungen für Werte und Ableitungen der Funktionen anzugeben. >>

Ich habe auch versucht, ihm zusätzliche Anfangsdaten zu geben, wie in der Fehlermeldung vorgeschlagen, aber ohne Erfolg. Die Frage Ich weiß nicht, ob es Potenzial gibt, meinen Code zu verbessern, oder ob ein Upgrade auf eine neuere Version von Mathematica das Problem lösen würde oder im schlimmsten Fall ein "zu hässliches" System für die numerische Behandlung ist.

Antworten

1 AlexTrounev Aug 26 2020 at 21:54

Um diese Art von Problemen zu lösen, können wir die Wellenfunktion in zwei Teile teilen $\psi=\psi_1+i\psi_2$. Wir verwenden auch einige Optionen NDSolve, um dieses Problem lösbar zu machen. Nehmen wir das an$\sigma$ ist echt, dann haben wir