# Nonlinear Schrodinger equation involving derivative of Abs causes ndnum warning

So, I'm using NDSolve to approximate the solution to a Schrödinger like equation

$i\psi_{t} + \psi_{xx} + \frac{1}{2}|\psi|^{2} \psi + \lambda \left( \left[|\psi|^{2} \psi\right]_{xx} + |\psi|^{4} + \psi \right) =0$

where $\psi: \mathbb{R}^{1+1} \to \mathbb{C}$. To do this I've modified the code from http://reference.wolfram.com/legacy(search for soliton to find the code) to include the $\lambda$ terms. I've tested the code for $\lambda =0$ and things are working just fine.

The problem I'm having is that when I include the second derivative term for $\lambda \neq 0$, Mathematica encounters a non-numerical derivative at $t=0$. I've looked at the second derivative of the initial condition

F[s_] = Soliton[1, π/12, m, 0][s, 0.]


and found that Mathematica yields

N[F[0]] = 2. ((0. +
0. I) + (0.999902 - 4.31242*10^-21 I) ((0. + 0.261799 I) -
0.261799 Derivative[1][Im][0.]))


and so I think that the problem is that it is keeping around the derivative of the imaginary part of a manifestly real quantity.

I'd like to get it to not work so hard at something so simple and would like to pass it something like $s \in \mathbb{R}$, but I can't figure out a way to do so. Is there a way to tell Mathematica to treat $s$ a little simpler?

Edit: Relevant codes.

Soliton[a_, b_, m_, s0_][s_, t_] := Module[{omega, ek, xx = s - s0},
omega = (2 m - 1) a^2 - b^2; Sqrt[m] a JacobiCN[a (xx - 2 b t), m] Exp[I (b xx + omega t)]];

per[m_] := 4 EllipticK[m];

m = mm /. FindRoot[per[mm] == 24, {mm, 0.5}];

ψ = u /.
First@ NDSolve[{I D[u[s, t], t] + D[u[s, t], {s, 2}] +
1/2 Abs[u[s, t]]^2  u[s, t] +
(1/16) (D[Abs[u[s, t]]^2 u[s, t], {s, 2}] +
Abs[u[s, t]]^4 u[s, t]) == 0,
u[s, 0.] == 2 Soliton[1, π/12, m, 0][s, 0],
u[-24, t] == u[24, t]}, u, {t, 0., 20.}, {s, -24., 24.},
MaxSteps -> Infinity, MaxStepSize -> 0.2, PrecisionGoal -> 5,
AccuracyGoal -> 5,
Method -> {"MethodOfLines", Method -> "StiffnessSwitching",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 100, "DifferenceOrder" -> "Pseudospectral"},
Method -> {"Adams", "MaxDifferenceOrder" -> 4}}];


NDSolve::ndnum

• I suspect Im'[..] comes from how you coded in Soliton and has little to do with s. You'd need to post the code for someone else to be able to check it out. May 23, 2017 at 20:03
• That's interesting. I've posted the codes and will next look into the issues associated with the module. Thanks! May 23, 2017 at 21:00
• Is m supposed to have a value? For N[F[0]] I get (1. + 0. I) Sqrt[m] JacobiCN[0., m] with what it posted. May 23, 2017 at 21:48
• Thanks for the update. Now I get 0.999951 - 4.31264*10^-21 I for N[F[0]] -- no Im'[]. Have you tried starting from a fresh kernel? May 24, 2017 at 2:02
• Yeah. I've included a code block that can be run straight-up. If you set the (1/16) term to 0 then you get a tidy traveling wave solution to the evolution. Keeping the 1/16 term throws the error, which goes away when you omit the second derivative term. May 24, 2017 at 2:12

The root of evil is Abs'[…], which is generated by D[Abs[u[s, t]]^2 u[s, t], {s, 2}] in the equation. (D cannot handle Abs[…] properly. ) To resolve the problem we need to express $|\psi|^{2}$ in a differentiable way. Notice that

$$\left| \psi \right| ^2=\psi \psi ^*$$ $$\frac{\partial \psi ^*}{\partial s}=\left(\frac{\partial \psi }{\partial s}\right)^*$$

We reform the equation to

eq = With[{abs = Abs[u[s, t]], abssquare = u[s, t] uconj[s, t]},
{I D[u[s, t], t] + D[u[s, t], {s, 2}] +
1/2 abs^2 u[s, t] + (1/16) (D[abssquare u[s, t], {s, 2}] + abs^4 u[s, t]) == 0,
u[s, 0.] == 2 Soliton[1, π/12, m, 0][s, 0],
u[-24, t] == u[24, t]}] /.
(uconj | Derivative[d__][uconj])[s, t] :> Conjugate[Derivative[d][u][s, t]];


You may find part of the definition of eq confusing. To understand it, just notice the following truth:

1. a /. a | b[m_] :> {m} outputs {}.
2. Derivative[][u] outputs u.

Now NDSolve can handle the equation with a manual adjustion for Method option:

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

sol = u /. First@
NDSolve[eq, u, {t, 0., 20}, {s, -24., 24.}, Method -> mol[51]]; // AbsoluteTiming
DensityPlot[Abs@sol[s, t], {s, -24, 24}, {t, 0, 20}]


At least in v9, the manual adjustion for Method is necessary, or NDSolve will spit out ndnum warning and fails again. I think it's another bug of NDSolve.