Coupled PDE system in atomic physics

Question

My question is about the implementation of a system of coupled PDEs to Mathematicas NDSolve routine. I'm considering one-dimensional toy model in atomic physics. The model describes two fields $\psi =\psi(t,z)$ and $\sigma= \sigma(z;t)$ coupled to each other i.e. $$ 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 $$ Additional I'm imposing periodic boundary conditions for $\psi(-L/2,t) = \psi(L/2,t)$ and $\sigma(-L/2,t) = \sigma(L/2,t)$ and set some initial conditions $\psi(z,0)=f(z)$ and $\sigma(z,0)=g(z)$.

EDITED:

Here is my current version of the 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}];

Unfortunately it comes with two problems, the first one concerns the Solver itself, since there is no time derivative in my equation for the second field $\sigma$ it handles the system as a DAE and give this two Warnings

NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations. >>

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. >>

I don't know if this is a "real" problem (I'm using Mathematica 9.x). The second one is more problematic, it concerns the amount of grid points used. This mainly comes from the equations itself I guess and cause an error that he can not find an appropriate solution within the tolerance bounds.

NDSolve::mxsst: Using maximum number of grid points 500 allowed by the MaxPoints or MinStepSize options for independent variable z. >>

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. >>

I also tried to give him additional initial data as suggested by the error message but without success. The question The thing I don't know is if there is any potential to improve my code, or if an upgrade to a newer version of Mathematica would solve the problem or in worst case its a "too ugly" system for numerical treatment.

Answer 1

To solve this kind of problems we can divide wave function into two parts $\psi=\psi_1+i\psi_2$. Also we use some options for NDSolve to make this problem solvable. Let suppose that $\sigma$ is real, then we have

(*constants*)h = 1;(*Planck constant*)m = 1;(*particle mass*)
Subscript[\[Alpha], 
  s] = 1;(*scattering length*)\[Omega] = 1;(*radial frequency*)V = \
0;(*longitudinal potential*)(*ranges*)L = 2;(*length of the box*)tmin \
= 0;
tmax = 0.1;

(*equations*)
eqn1 = { D[\[Psi]1[z, t], t] == -h^2/(2 m) D[\[Psi]2[z, t], z, z] + 
     V \[Psi]2[z, t] + 
     h^2 Subscript[\[Alpha], s]/
       m \[Sigma][z, 
        t]^(-2) (\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2) \[Psi]2[z, t] + 
     h^2/(2 m) \[Sigma][z, t]^(-2) \[Psi]2[z, t] + 
     m \[Omega]/2 \[Sigma][z, t]^2 \[Psi]2[z, t] + 
     h^2/(2 m) \[Sigma][z, t]^(-2) D[\[Sigma][z, t], z]^2 \[Psi]2[z, 
       t], - D[\[Psi]2[z, t], 
      t] == -h^2/(2 m) D[\[Psi]1[z, t], z, z] + V \[Psi]1[z, t] + 
     h^2 Subscript[\[Alpha], s]/
       m \[Sigma][z, 
        t]^(-2) (\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2) \[Psi]1[z, t] + 
     h^2/(2 m) \[Sigma][z, t]^(-2) \[Psi]1[z, t] + 
     m \[Omega]/2 \[Sigma][z, t]^2 \[Psi]1[z, t] + 
     h^2/(2 m) \[Sigma][z, t]^(-2) D[\[Sigma][z, t], z]^2 \[Psi]1[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]/(\[Psi]1[z, t]^2 + \[Psi]2[z, 
          t]^2) (D[(\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2), 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) (\[Psi]1[z, t]^2 + \[Psi]2[z, t]^2);

(*boundary conditions*)
bc = {\[Psi]1[L/2, t] == \[Psi]1[-L/2, t], \[Psi]2[L/2, 
     t] == \[Psi]2[-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]1[z, 0] == Exp[-((z)^2)], \[Psi]2[z, 0] == 0};
(*solve system*)
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{Psi1, Psi2, S} = 
    NDSolveValue[{eqn1, eqn2, bc, bcwidth, icwave, 
      icwidth}, {\[Psi]1, \[Psi]2, \[Sigma]}, {z, -L/2, L/2}, {t, 
      tmin, tmax}, 
     Method -> {"IndexReduction" -> Automatic, 
       "EquationSimplification" -> "Residual", 
       "PDEDiscretization" -> {"MethodOfLines", 
         "SpatialDiscretization" -> {"TensorProductGrid", 
           "MinPoints" -> 81, "MaxPoints" -> 81, 
           "DifferenceOrder" -> "Pseudospectral"}}}, 
     EvaluationMonitor :> (currentTime = t;)];]; 

Visualization of numerical solution

{Plot3D[Psi1[z, t], {z, -L/2, L/2}, {t, tmin, tmax}, Mesh -> None, 
  ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "Re\[Psi]"], 
 Plot3D[Psi2[z, t], {z, -L/2, L/2}, {t, tmin, tmax}, Mesh -> None, 
  ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "Im\[Psi]"], 
 Plot3D[S[z, t], {z, -L/2, L/2}, {t, tmin, tmax}, Mesh -> None, 
  ColorFunction -> Hue, AxesLabel -> Automatic, 
  PlotLabel -> "\[Sigma]", PlotRange -> All]}