# MoL: How to enforce Chebyshev–Gauss–Lobatto points in SpatialDiscretization?

In Mathematica documentation one is prompted to use a grid with points at the zeros of the Chebyshev polynomials so that Runge's phenomena arising from DifferenceOrder->"Pseudospectral" are eliminated.

Although there is an explicit example of how to construct such a grid

 CGLGrid[x0_, L_, n_Integer /; n > 1] := x0 + 1/2 L (1 - Cos[π Range[0, n - 1]/(n - 1)])
cgrid = CGLGrid[-5, 10.,16];


I have not yet found how to introduce this grid in SpatialDiscretization option.

Can anyone post a simple example?

I'm afraid you've misunderstood the document. The document actually means, when DifferenceOrder->"Pseudospectral" is chosen for non-periodic b.c., Chebyshev–Gauss–Lobatto (CGL) grid will be automatically used so that Runge's phenomena won't be extreme. This can be verified by

points = 35;

usol = NDSolveValue[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0, u[t, 0] == Sin[t],
u[t, 5] == 0}, u, {t, 0, 10}, {x, 0, 5},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> points,
"MinPoints" -> points, "DifferenceOrder" -> "Pseudospectral"}}];

xcoord = usol["Coordinates"][[2]];

CGLGrid[x0_, L_, n_Integer /; n > 1] := x0 + 1/2 L (1 - Cos[Pi Range[0, n - 1]/(n - 1)])

cgrid = CGLGrid[0, 5., points];

xcoord == cgrid
(* True *)


Still, you can use CGL grid for other difference order as shown in user21's answer, but I doubt if CGL grid helps in those cases. (If CGL grid is really a catholicon, why isn't it the default setting of NDSolve? )

• Nothing to say Mathematica wise, but I was amused by the word "catholicon". Can’t say I’ve ever seen or heard that word, even in some pretty old texts. Appears to just mean panacea, which is the common word, but always fun to learn new ones. – b3m2a1 Sep 7 at 20:10
• @b3m2a1 Oh, didn't know this word is deserted. I just checked the dictionary and copy and paste it :D . – xzczd Sep 8 at 5:12
• I figured as much :) archaic words like that mostly pop up in translation, esp. from languages with some distance from the target language. Still fun to learn. – b3m2a1 Sep 8 at 5:13
• I think there are two issues, choosing the grid and choosing a method of approximating derivatives. With "DifferenceOrder" -> n, for an integer n, one is effectively using piecewise interpolation over n+1 point subgrids for the "finite difference derivatives." The use of a Chebyshev grid will hardly make a difference, unless n is the full order. The "Pseudospectral" setting uses the full order plus the FFT to compute derivatives (the docs imply), which will perform better than the dense finite-difference-derivative matrices. – Michael E2 Oct 6 at 12:33
• @MichaelE2 Er…do you mean that "SpatialDiscretization" -> {"TensorProductGrid", "Coordinates" -> cgrid, "DifferenceOrder" -> points} is mathematically equivalent to "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> points, "MinPoints" -> points, "DifferenceOrder" -> "Pseudospectral", except that the algorithm used in the former case is less efficient ? – xzczd Oct 7 at 5:12

If you read the linked tutorial to the end you will find the following code:

mygrid =  Join[-5. + 10 Range[0, 48]/80, 1. + Range[1, 4 70]/70];
ν = 0.01;
bsolg = First[
NDSolve[{D[u[x, t], t] == ν D[u[x, t], x, x] -
u[x, t] D[u[x, t], x], u[x, 0] == E^-x^2, u[-5, t] == u[5, t]},
u, {x, -5, 5}, {t, 0, 4},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"Coordinates" -> {mygrid}}}]]