Système PDE couplé en physique atomique

Aug 19 2020

Ma question concerne la mise en œuvre d'un système de PDE couplés à la routine Mathematicas NDSolve. J'envisage un modèle de jouet unidimensionnel en physique atomique. Le modèle décrit deux champs$\psi =\psi(t,z)$ et $\sigma= \sigma(z;t)$ couplés les uns aux autres ie $$ 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 $$ J'impose des conditions aux limites périodiques supplémentaires pour $\psi(-L/2,t) = \psi(L/2,t)$ et $\sigma(-L/2,t) = \sigma(L/2,t)$ et fixez quelques conditions initiales $\psi(z,0)=f(z)$ et $\sigma(z,0)=g(z)$.

ÉDITÉ:

Voici ma version actuelle du code

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

Malheureusement, cela pose deux problèmes, le premier concerne le solveur lui-même, car il n'y a pas de dérivée temporelle dans mon équation pour le deuxième champ $\sigma$ il gère le système comme un DAE et donne ces deux avertissements

NDSolve :: pdord: Certaines des fonctions ont un ordre différentiel nul, donc les équations seront résolues comme un système d'équations différentielles-algébriques. >>

NDSolve :: mconly: Pour la méthode IDA, seul le code réel de la machine est disponible. Impossible de continuer avec des valeurs complexes ou au-delà des exceptions à virgule flottante. >>

Je ne sais pas si c'est un "vrai" problème (j'utilise Mathematica 9.x). Le second est plus problématique, il concerne la quantité de points de grille utilisés. Cela vient principalement des équations elles-mêmes je suppose et provoque une erreur selon laquelle il ne peut pas trouver une solution appropriée dans les limites de tolérance.

NDSolve :: mxsst: Utilisation du nombre maximal de points de grille 500 autorisé par les options MaxPoints ou MinStepSize pour la variable indépendante z. >>

NDSolve :: icfail: impossible de trouver les conditions initiales qui satisfont la fonction résiduelle dans les tolérances spécifiées. Essayez de donner des conditions initiales pour les valeurs et les dérivées des fonctions. >>

J'ai également essayé de lui donner des données initiales supplémentaires comme suggéré par le message d'erreur mais sans succès. La question La chose que je ne sais pas, c'est s'il y a un potentiel pour améliorer mon code, ou si une mise à niveau vers une version plus récente de Mathematica résoudrait le problème ou dans le pire des cas, c'est un système "trop laid" pour un traitement numérique.

Réponses

1 AlexTrounev Aug 26 2020 at 21:54

Pour résoudre ce genre de problèmes, nous pouvons diviser la fonction d'onde en deux parties $\psi=\psi_1+i\psi_2$. Nous utilisons également certaines options pour NDSolverésoudre ce problème. Supposons que$\sigma$ est réel, alors nous avons