2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ "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? $\endgroup$ – J. M.'s technical difficulties Nov 3 '17 at 3:03
  • $\begingroup$ 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. $\endgroup$ – Bill Watts Nov 3 '17 at 4:19
2
$\begingroup$

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]

enter image description here

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]

enter image description here

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]

enter image description here

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. enter image description here

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.

| improve this answer | |
$\endgroup$
0
$\begingroup$

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

enter image description here

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

enter image description here

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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}$ $\endgroup$ – Buddhapus Nov 4 '17 at 6:40
  • $\begingroup$ @Buddhapus You were right about $x^{-1}$ but I don't think anything else is wrong. See the edited version. $\endgroup$ – MathX Nov 6 '17 at 17:23

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.