Sistema PDE accoppiato in fisica atomica

Aug 19 2020

La mia domanda riguarda l'implementazione di un sistema di PDE accoppiati alla routine Mathematicas NDSolve. Sto considerando un modello giocattolo unidimensionale in fisica atomica. Il modello descrive due campi$\psi =\psi(t,z)$ e $\sigma= \sigma(z;t)$ accoppiati tra loro, ad es $$ 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 $$ Inoltre sto imponendo condizioni al contorno periodiche per $\psi(-L/2,t) = \psi(L/2,t)$ e $\sigma(-L/2,t) = \sigma(L/2,t)$ e impostare alcune condizioni iniziali $\psi(z,0)=f(z)$ e $\sigma(z,0)=g(z)$.

MODIFICATO:

Ecco la mia versione attuale del codice

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

Sfortunatamente presenta due problemi, il primo riguarda il Risolutore stesso, poiché non esiste una derivata temporale nella mia equazione per il secondo campo $\sigma$ gestisce il sistema come DAE e fornisce questi due avvisi

NDSolve :: pdord: Alcune delle funzioni hanno ordine differenziale zero, quindi le equazioni verranno risolte come un sistema di equazioni algebriche differenziali. >>

NDSolve :: mconly: per il metodo IDA, è disponibile solo il codice reale della macchina. Impossibile continuare con valori complessi o oltre le eccezioni a virgola mobile. >>

Non so se questo sia un problema "reale" (sto usando Mathematica 9.x). Il secondo è più problematico, riguarda la quantità di punti della griglia utilizzati. Questo deriva principalmente dalle equazioni stesse immagino e causa un errore che non riesce a trovare una soluzione adeguata entro i limiti di tolleranza.

NDSolve :: mxsst: utilizzo del numero massimo di punti della griglia 500 consentito dalle opzioni MaxPoints o MinStepSize per la variabile indipendente z. >>

NDSolve :: icfail: impossibile trovare le condizioni iniziali che soddisfano la funzione residua entro le tolleranze specificate. Prova a fornire le condizioni iniziali sia per i valori che per le derivate delle funzioni. >>

Ho anche provato a fornirgli ulteriori dati iniziali come suggerito dal messaggio di errore ma senza successo. La domanda La cosa che non so è se c'è un potenziale per migliorare il mio codice, o se un aggiornamento a una versione più recente di Mathematica risolverebbe il problema o nel peggiore dei casi è un sistema "troppo brutto" per il trattamento numerico.

Risposte

1 AlexTrounev Aug 26 2020 at 21:54

Per risolvere questo tipo di problemi possiamo dividere la funzione d'onda in due parti $\psi=\psi_1+i\psi_2$. Inoltre utilizziamo alcune opzioni per NDSolverendere risolvibile questo problema. Supponiamo che$\sigma$ è reale, allora abbiamo