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.
icdw = D[σ[z], z][0] == 0;
is obviously wrong, if you mean $σ'(0)=0$, it should beicdw = σ'[0] == 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 $σ$? $\endgroup$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 beicdwidth = D[\[Sigma][z, t], z] == 0 /. z -> 0;
. But I believeicdwidth
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. $\endgroup$