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}}];


  • 2
    $\begingroup$ 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. $\endgroup$ – Michael E2 May 23 '17 at 20:03
  • $\begingroup$ That's interesting. I've posted the codes and will next look into the issues associated with the module. Thanks! $\endgroup$ – user153764 May 23 '17 at 21:00
  • $\begingroup$ 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. $\endgroup$ – Michael E2 May 23 '17 at 21:48
  • $\begingroup$ 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? $\endgroup$ – Michael E2 May 24 '17 at 2:02
  • $\begingroup$ 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. $\endgroup$ – 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}]

Mathematica graphics

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.