NSolve Slow In Finding Solutions to Nonlinear Equations

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.

• You have 8 equations (Length[list]) and 7 variables.
– JimB
Apr 1, 2017 at 2:11
• With your edit did you change your code listed above at all?
– JimB
Apr 1, 2017 at 5:42
• @JimBaldwin Yes, just changed it! Should produce an empty set now. Apr 1, 2017 at 5:46
• Rather than NSolve how about FindRoot: sol = FindRoot[ list, {{c[1], 0}, {c[2], 0}, {c[3], 0}, {c[4], 0}, {c[5], 0}, {c[6], 0}, {c[7], 0}}]
– JimB
Apr 1, 2017 at 5:53
• No genius here. I've just made many more programming errors than you and have remembered at least a few of those to know what to look for.
– JimB
Apr 1, 2017 at 17:54

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 *)


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

• How to plot solution u[x] ? Mar 25, 2018 at 7:35
• @MariuszIwaniuk Something like chebF[c_, x_] := c.Cos[2 Range[0, Length@c - 1] ArcCos[Sqrt@x]]; usol = chebF[sols[[-1, ;; -20]], x]; Plot[usol, {x, 0, 1}]. Mar 25, 2018 at 10:23
• If γ0 = 0;β0 = 1;α = 2;fde = -sum0[m] + s[x, M] - 2 (s[x, M])^3 comparison with NDSolve[{-u[x] + u''[x] - 2*u[x]^3 == 0, u[0] == 0, u'[0] == 1}, u, {x, 0, 2}]` plots not the same? Mar 25, 2018 at 10:53
• Can your approach be used for FDE bvp? If yes then can you please give some hints to solve equations 6.1 and 6.2 in this link? sciencedirect.com/science/article/pii/S1110016816000156
– zhk
Oct 26, 2019 at 1:38