Sistema PDE acoplado en física atómica

Aug 19 2020

Mi pregunta es sobre la implementación de un sistema de PDE acoplados a la rutina NDSolve de Mathematicas. Estoy considerando un modelo de juguete unidimensional en física atómica. El modelo describe dos campos$\psi =\psi(t,z)$ y $\sigma= \sigma(z;t)$ acoplados entre sí, es decir $$ 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 $$ Adicional Estoy imponiendo condiciones de contorno periódicas para $\psi(-L/2,t) = \psi(L/2,t)$ y $\sigma(-L/2,t) = \sigma(L/2,t)$ y establece unas condiciones iniciales $\psi(z,0)=f(z)$ y $\sigma(z,0)=g(z)$.

EDITADO:

Aquí está mi versión actual del código

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

Desafortunadamente, viene con dos problemas, el primero concierne al Solver en sí mismo, ya que no hay derivada de tiempo en mi ecuación para el segundo campo. $\sigma$ maneja el sistema como un DAE y da estas dos advertencias

NDSolve :: pdord: Algunas de las funciones tienen orden diferencial cero, por lo que las ecuaciones se resolverán como un sistema de ecuaciones algebraicas diferenciales. >>

NDSolve :: mconly: Para el método IDA, solo está disponible el código real de la máquina. No se puede continuar con valores complejos o más allá de las excepciones de punto flotante. >>

No sé si esto es un problema "real" (estoy usando Mathematica 9.x). El segundo es más problemático, se refiere a la cantidad de puntos de cuadrícula utilizados. Esto proviene principalmente de las ecuaciones en sí, supongo, y provocan un error que no puede encontrar una solución adecuada dentro de los límites de tolerancia.

NDSolve :: mxsst: Usando el número máximo de puntos de cuadrícula 500 permitidos por las opciones MaxPoints o MinStepSize para la variable independiente z. >>

NDSolve :: icfail: No se pueden encontrar las condiciones iniciales que satisfacen la función residual dentro de las tolerancias especificadas. Intente dar condiciones iniciales tanto para los valores como para las derivadas de las funciones. >>

También traté de darle datos iniciales adicionales como lo sugiere el mensaje de error, pero sin éxito. La pregunta Lo que no sé es si hay algún potencial para mejorar mi código, o si una actualización a una versión más nueva de Mathematica resolvería el problema o, en el peor de los casos, es un sistema "demasiado feo" para el tratamiento numérico.

Respuestas

1 AlexTrounev Aug 26 2020 at 21:54

Para resolver este tipo de problemas podemos dividir la función de onda en dos partes $\psi=\psi_1+i\psi_2$. También usamos algunas opciones para NDSolveresolver este problema. Supongamos que$\sigma$ es real, entonces tenemos