NSolve Slow In Finding Solutions to Nonlinear Equations

Question

I'm working on numerically solving the fractional nonlinear Schrödinger equation (here of order $\alpha$):

$$\frac{1}{2}u-\frac{1}{2}\frac{\partial^{\alpha}}{\partial x^{\alpha}}u-u^3=0.$$

The paper I read uses shifted Chebyshev polynomials to construct a series solution to any fractional nonlinear differential equation [1]. They are shifted from their regular domain $[-1,1]$ to $[0,1]$. Essentially what we solve for are the coefficients in front of each shifted Chebyshev polynomial, here array c, with 7 total coefficients. I think I have the code set up right, but NSolve has run for about thirty minutes now and is still working. I'm afraid that it's not going to find the coefficients because of the nonlinearity. Is there anything I can do to improve the performance or ensure the coefficients are found? So far no errors. Perhaps there is a different way of solving for them.

ClearAll["Global`*"]
γ0 = 0; (* u(0) = 0 *)
β0 = 0.8; (* u'(0) = 0.8 *)
α = 18/15;
m = 6; (* number of terms in approximation series-1 *)

Array[c, m + 1]; (* array of coefficients *)

ts[n_, x_] := ChebyshevT[n, 2 x - 1]; (* Defines Shifted Chebyshev Polynoms *)

pts = N[Solve[{ts[5, x] == 0}, x]]
collocation = Table[pts[[i, 1, 2]], {i, 1, Length[pts]}] (* collocation points *)

s[x_] := Sum[c[n] ts[n - 1, x], {n, 1, m + 1}]; (* series approximation *)

b[n_, r_, α_] := (-1)^r 2^(2 n - 2 r) n (Factorial[2 n - r - 1] Factorial[n - r])/
(Factorial[r] Factorial[2 n - 2 r] Gamma[n - r + 1.0 - α]);

eqn1 = Sum[(-1)^(n - 1) c[n], {n, 1, m + 1}] == γ0; (*u(0) initial condition*)
eqn2 = 2 Sum[(-1)^n (n - 1)^2 c[n], {n, 1, m + 1}] == β0; (*u'(0) initial condition*)
eqn[x_] := Sum[Sum[c[n + 1] b[n, r, α] x^(n - r - α), {r, 0, n - 2}], {n, 2, m}] == 
s[x] - 2 (s[x])^3 (* fractional nonlinear Schroedinger *)
eqnlist = Table[eqn[x] /. x -> collocation[[i]], {i, 1, Length[collocation]}];
(* fractional nonlinear Schroedinger evaluated at collocation points *)

list = Flatten[{eqn1, eqn2, eqnlist}];

NSolve[list, {c[1], c[2], c[3], c[4], c[5], c[6], c[7]}, Reals]

Thanks for any help!

Edit Per comment below I changed the number of collocation points to 5 (used to be 6). The nonlinear Schrödinger equation evaluated at 5 collocation points, plus 2 equations for initial conditions, gives us 7 total equations, for 7 unknowns. Unfortunately Mathematica gives me an empty set.

---

[1]: A Chebyshev Pseudo-Spectral Method for Solving Fractional-Order Integro-Differential Equations

Answer 1

Changing from NSolve to FindRoot gives a solution:

sol = FindRoot[list, {{c[1], 0}, {c[2], 0}, {c[3], 0}, {c[4], 0}, {c[5], 0}, c[6], 0}, {c[7], 0}}]
{c[1] -> 0.45581, c[2] -> 0.459006, c[3] -> -0.0111779, c[4] -> -0.0173396, c[5] -> -0.00224159, 
 c[6] -> 0.000872271, c[7] -> 0.000148173}

As a check on the solution:

list /. sol
(* {False,True,True,True,True,True,True} *)

We see that the first element list is False. But further examination of that first element shows that isn't a problem. The first element is

list[[1]]
(* c[1] - c[2] + c[3] - c[4] + c[5] - c[6] + c[7] == 0 *)
c[1] - c[2] + c[3] - c[4] + c[5] - c[6] + c[7] /. sol
(* 9.54098*10^-18 *)

Answer 2

One thing to watch out for is that the formulation of the equations is not particularly numerically stable. The formula for eqn[x] in particular has large coefficients that lead to substantial precision loss, presumably due to subtraction of like sized terms.

The precision loss can be controlled somewhat by using exact coefficients (e.g. β = 8/10) and by evaluating b[n, r, α] and one of the sums in eqn[x] symbolically and simplifying them:

Block[{n, r, α},
 b[n_, r_, α_] =
  (-1)^r 2^(2 n - 2 r) n (Factorial[2 n - r - 1] Factorial[n - r]) / 
   (Factorial[r] Factorial[2 n - 2 r] Gamma[n - r + 1 - α]) // FullSimplify
 ]
Block[{m, n},
 sum0[m_] = Sum[Sum[c[n + 1] b[n, r, α] x^(n - r - α), {r, 0, n - 2}], {n, 2, m}]]

There may be other approaches, too. These things, however only get you so far. You can get almost single precision up to order 15-17 (see below). To get more, arbitrary precision is needed.

In the following, sys[M, initWP] constructs the OP's system of order M with an "initial working precision" initWP, which is the precision at which the Chebyshev points are calculated (directly using trigonometry, instead of solving for the roots of the polynomials). I use FindRoot like @Jim Baldwin (NSolve does work on the OP's problem if exact coefficients are used, but FindRoot is much faster.) The solver solve takes a vector of initial guesses for the Chebyshev coefficients; the length of the vector determines the order of the approximation sought.

ClearAll[ts, s, b, sys, c, solve]
SetAttributes[c, NHoldAll];
γ0 = 0;(*u(0)=0*)
β0 = 8/10;(*u'(0)=0.8*)
α = 18/15; 

ts[n_, x_] := ChebyshevT[n, 2 x - 1];(*Defines Shifted Chebyshev Polynoms*)
s[x_, M_] := Sum[c[n] ts[n - 1, x], {n, 1, M}];(*series approximation*)
Block[{n, r, α},
 b[n_, r_, α_] =
  (-1)^r 2^(2 n - 2 r) n (Factorial[2 n - r - 1] Factorial[n - r]) / 
   (Factorial[r] Factorial[2 n - 2 r] Gamma[n - r + 1 - α]) // FullSimplify
 ];
Block[{m, n},
 sum0[m_] = Sum[Sum[c[n + 1] b[n, r, α] x^(n - r - α), {r, 0, n - 2}], {n, 2, m}]];
sys[M_Integer, initWP_: 50] := 
  Module[{m = M - 1, n = M - 2, collocation, ics, fde},
   collocation = 
    N[Rescale[Sin[Range[-(n - 1)/2, (n - 1)/2] Pi/n], {-1, 1}, {0, 1}], initWP];
   ics = {Sum[(-1)^(n - 1) c[n], {n, 1, M}] == γ0,  (*u(0) initial condition*)
     2 Sum[(-1)^n (n - 1)^2 c[n], {n, 1, M}] == β0  (*u'(0) initial condition*)};
   fde = -sum0[m] + s[x, M] - 2 (s[x, M])^3;   (*fractional nonlinear Schroedinger (FNLS)*)
   fde = Thread[(fde /. x -> collocation) == 0] /. 
     z_Real /; z == 0 :> 0;                    (*FNLS at collocation pts*)

   Flatten[{ics, fde}]
   ];

solve[c0_List, initWP_: 50] :=  
  With[{M = Length@c0, vars = Array[c, Length@c0]},
   With[{vars0 = Transpose[{vars, c0}]},
    vars /. FindRoot[sys[M, initWP], vars0]
    ]];

OP's example:

solve[ConstantArray[0, 5]]
(*  {0.457648, 0.462563, -0.00862945, -0.017091, -0.00354613}  *)

If things converge nicely, then for the computed Chebyshev expansion and the Cheybshev series of the true solution, respectively $$\tilde u(x) = \sum_{j=1}^M c_j T_{j-1}(x), \quad u(x) = \sum_{j=1}^\infty a_j T_{j-1}(x)\,,$$ we should have $c_j \rightarrow a_j$ as $M -> \infty$. Further for differentiable function, $a_j \rightarrow 0$ rapidly. Thus we can get a picture of the quality of the approximation, by plotting the Chebyshev coefficients.

Here is the plot of the solutions up to M == 102 (this takes a minute or so to compute):

PrintTemporary@Dynamic@foo;
workprec = 100;
sols = Rest@NestList[
    solve[foo = Length@#; Join[#, {0}], workprec] &,
    ConstantArray[0, 3], 100];

ListLinePlot[RealExponent@sols, 
 PlotStyle -> 
  Table[Directive[Opacity[0.5], ColorData["Rainbow"][i/Length@sols]], {i, Length@sols}],
 PlotLabel -> Row[{"Precision: sys -> ", Precision[sys[Length@sols, workprec]]}],
 Frame -> True, FrameLabel -> {"order n", HoldForm@Log10@Abs@Subscript[c, n]}
 ]

There is a divergence at the tail that is due to the nonlinearity. If one drops the last 20 or so coefficients, one appears to get a series that is a pretty good approximation. The Chebyshev series is converging somewhat slowly (for a Chebyshev series), so getting double precision will be difficult, but probably feasible. One can also see that in the construction of the last system, the precision loss is about 85 digits! Roughly speaking one wants the initial working precision to be about 20 + 3/4 * M (the 3/4 seems a little too high).

By way of comparison, below on the left is the computation with workprec = 20, and on the right, is the OP's set up repeated for different orders. The OP's converges to the high precision solution for M up to about 25 and calculates the coefficients up to order 15 or so accurately. After that, it diverges