# 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. – Michael E2 May 23 '17 at 20:03
• That's interesting. I've posted the codes and will next look into the issues associated with the module. Thanks! – user153764 May 23 '17 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. – Michael E2 May 23 '17 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? – Michael E2 May 24 '17 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. – user153764 May 24 '17 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.