# Coupled PDE system in atomic physics

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.

• Have you read this?: mathematica.stackexchange.com/q/188197/1871 Aug 20, 2020 at 5:55
• Yes, but the problem is, the most of the questions are unanswered and among those are answered, I couldn't find a suitable solution for me. Aug 20, 2020 at 7:50
• At least 3 issues: 1. icdw = D[σ[z], z] == 0; is obviously wrong, if you mean $σ'(0)=0$, it should be icdw = σ' == 0;, 2. If $σ$ varies over time, then this is not a coupled system of PDE and ODE, but a PDAE system, and σ[z] should be σ[z, t]. We already have some related posts e.g. this: mathematica.stackexchange.com/q/184281/1871 3. What's the b.c. for $σ$? Aug 22, 2020 at 15:42
• @xzczd Hi, I followed your advice and changed the dependence of σ (see the edited code). Of course, you are right and it should of course depend on time that NDSolve can solve it in every timestep, and added the corresponding bc. Unfortunately, it didn't change the error, maybe this could be a discretization problem as well, as you pointed out in your answer on the given link above? Aug 24, 2020 at 14:12
• Still, icdwidth = D[\[Sigma][z, t], z][0, t] == 0; is obviously wrong. If you mean $\frac{\partial \sigma}{\partial z}\big|_{z=0}=0$, then it should be icdwidth = D[\[Sigma][z, t], z] == 0 /. z -> 0;. But I believe icdwidth is just redundant, because you've already imposed periodic b.c. in $z$ direction. Even in $t$ direction, you just need at least one i.c. for $\sigma$, because there's no derivative of $\sigma$ in $t$ direction in the system. Aug 25, 2020 at 5:08

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]} • Hi, your answer is almost perfect, the plots look very promising. But I would like to ask how long the code took to run? Because I'm running it on the university cluster and it seems to stuck somewhere, but it isn't producing any error messages. Could this be a version conflict or something which includes additional updates for Mathematica? Aug 27, 2020 at 12:04
• @Hamilcar What version do you run? My version is 12+ for Windows 10. But actually this code can be use on version 11+ too. Aug 27, 2020 at 12:12
• I'm running version 9.0 Linux/Windows, but I already contacted the IT department for a version update. Aug 27, 2020 at 12:52
• @Hamilcar With number of points 81 code takes 9.73708 second on my computer. May be better you run just NDSolve without Dynamic["time: " <> ToString[CForm[currentTime]]] and AbsoluteTiming[]. Aug 27, 2020 at 14:37
• Hi, the code worked in the new version at a similar time. The problem is solved Aug 27, 2020 at 15:08