# Fractional nonlinear Schrodinger equation … looking for cnoidal solution

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.

• "Since the first derivative should be undefined..." - NDSolve[] will definitely have trouble. Recall that it uses piecewise polynomials for interpolation, and that polynomials have never had a vertical slope. Consider factoring out the singular behavior so that you have a solution that can be approximated by piecewise polynomials, and then multiply that isolated factor with the resulting interpolant? – J. M.'s torpor Nov 3 '17 at 3:03
• Another problem I see is that if you take your pde, multiply it by sqrt[x] to get rid of the singularity at x=0, and then set x to 0 in the resulting pde, you get u[0] = 0. So the pde wants u[0] = 0, and you are forcing u[a] to be 1 as a approaches 0, making instability expected. – Bill Watts Nov 3 '17 at 4:19

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.

I think any numerical differentiation model will have an issue with your singularity and blowup at $x=0$ unless you reformulate your problem somehow. For your simplified model where you are sure $x>0$, you can multiply both side of your differential equation in $x^{q-2}$ and solve that instead. Your original equation is:

$-\frac {(q - 1) q x^{2 - q} u'' (x)} {2\Gamma (3 - q)} - \frac {q x^{1 - q} u' (x)} {\Gamma (2 - q)} + u ( x)\left (1 - \frac {x^{-q}} {\Gamma (1 - q)} \right) - 2 u \left( x\right)^3=0$

After multiplying each term by $x^{q-2}$:

$-\frac {(q - 1) q u'' (x)} {2\Gamma (3 - q)} - \frac {q u' (x)} {x \Gamma (2 - q)} + u ( x)\left (x^{q - 2} - \frac {x^{2}} {\Gamma (1 - q)} \right) - 2 x^{q-2} u \left( x\right)^3=0$

sol[q_, a_] :=
NDSolveValue[{(x^(-2 + q) - x^2/Gamma[1 - q]) u[x] -
2 x^(-2 + q) u[x]^3 - (q u'[x])/(
x Gamma[2 - q]) - ((-1 + q) q u''[x])/(2 Gamma[3 - q]) == 0,
u[a] == 1, u'[a] == 0}, u[x], {x, a, 100}]


For your values of $q$ and $a$ the solution near zero is

If $a=10^{-7}$ we can still see the same behavior:

Seems to me that you will always get the decaying solution you want near zero. However, I think you are basically changing the mathematics of your problem by considering $x>0$ and it is not what you might want in reality so this simplification trick could be the issue.

• I think you would have $x^{-2}$ instead of $x^2$ and $x^{-1}$ instead of $x$, if you're multiplying everything by $x^{q-2}$ – Buddhapus Nov 4 '17 at 6:40
• @Buddhapus You were right about $x^{-1}$ but I don't think anything else is wrong. See the edited version. – MathX Nov 6 '17 at 17:23