Fractional nonlinear Schrodinger equation ... looking for cnoidal solution

Question

I need to solve the following gnarly differential equation (a version of fractional nonlinear Schrödinger equation):

q = 3/2.;
a = 0.0001;
sol = NDSolveValue[{(1 - x^(-q)/Gamma[1 - q]) u[x] - (q x^(1 - q))/Gamma[2 - 
q] u'[x] - (q (q - 1) x^(2 - q))/(2 Gamma[3 - q]) u''[x] - 2 (u[x])^2 u[x] == 0,
u[a] == 1,
u'[a] == 0},
u[x], {x, a, 100}]

Due to there being $x^{-q}$ in the equation, $x > 0$. However, my initial condition is defined in terms of $x = 0$ so the way I fix it is by closely approaching the origin from the right at $x = a$.

Then I plot the solution and look at it. Here's the problem: the solution changes as I make a smaller and smaller. It makes me worried that whatever Mathematica is getting is not convergent, but perhaps the problem is not well-posed either.

The reason is that the solution is supposed to be like a cnoidal wave (whose shape is a function of q), which develops an undefined first slope u'[x] at the origin. Since the first derivative should be undefined, then it is debilitating to specify it in the NDSolve command. The solution is also supposed to go to 1 at $x = 0$ but as I vary a, the solution seems to diverge in the region $a \leq x \leq 0.02$. From my basic understanding of the problem, the solution should start at 1 and decay from there.

Has anyone ran into this problem before (maybe in fluid simulations)? I'm not sure how to address these problems in Mathematica.

Answer 1

It looks like you really need to use a small a to see what is going on with this system.

pde = (1 - x^(-q)/Gamma[1 - q]) u[x] - (q x^(1 - q))/Gamma[2 - q] u'[x] -
(q (q - 1) x^(2 - q))/(2 Gamma[3 - q]) u''[x] - 
       2 (u[x])^2 u[x] == 0;

Use your values for u and a fairly small a.

q = 3/2
a = 10^-11
Clear[u]

Increase WorkingPrecision and MaxSteps such that

sol = NDSolve[{pde2, u[a] == 1, u'[a] == 0}, u[x], {x, a, 200}, 
  WorkingPrecision -> 40, MaxSteps -> 500000];

u[x_] = u[x] /. sol[[1]]

strg = {"a=", N[a]}[[1]] + {"a=", N[a]}[[2]];
Plot[u[x], {x, a, 200}, PlotRange -> {-1, 1}, PlotLabel -> strg]

Its a little hard to see, but for small x, u oscillates about u = 0 and then jogs to oscillating about a positive u. I won't show it here, but if I change to a = 10^-12, the oscillations jog to a negative u. Check when a becomes really small:

q = 3/2
a = 10^-22
Clear[u]

sol = NDSolve[{pde2, u[a] == 1, u'[a] == 0}, u[x], {x, a, 200}, 
  WorkingPrecision -> 40, MaxSteps -> 500000];

u[x_] = u[x] /. sol[[1]]

strg = {"a=", N[a]}[[1]] + {"a=", N[a]}[[2]];
Plot[u[x], {x, a, 200}, PlotRange -> {-10, 10}, PlotLabel -> strg]

In this case, the oscillations remain about u = 0 until around x = 150 This may be the behavior you are after. Smaller values of a are possible, but NDSolve starts failing for this system when WorkingPrecision gets too high. Looking at the behavior near 0:

Plot[u[x], {x, a, 2 a}, PlotRange -> All, PlotLabel -> strg]

This shows that of u'[0]==0 at x = a, at least by eyeball. The next plot shows why getting rid of the singularity by multiplying the pde by x^3/2 does not help.

This shows a huge peak at small x, but still a positive x. The closer to zero a gets, the higher the peak, which is evidently why NDSolve reports an infinity for the equation with the singularity removed and using a = 0. It is also true that if you remove the singularity from the pde, solving the subsequent pde for u''[x] still has x^2 in the denominator.