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

Question

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?

score 5 · Accepted Answer · answered Sep 06 '19 at 11:36

5

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

answered Sep 06 '19 at 11:36

xzczd

65,995
9
163
468

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 07 '19 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 08 '19 at 05: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 08 '19 at 05: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 06 '19 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 07 '19 at 05:12
Yes, I think that's correct, though perhaps it's "DifferenceOrder" -> points - 1. (But numerically and with respect to speed, "Pseudospectral" should perform better.) – Michael E2 Oct 07 '19 at 09:46
@M The correct DifferenceOrder for FiniteDifferenceDerivative seems to be points-d, where d is derivative order. (FiniteDifferenceDerivative automatically reduces the order when DifferenceOrder is too high, though. ) Test: grid = N[CGLGrid[0,5,points],32];d={2,3};dm=NDSolve\FiniteDifferenceDerivative[d, {grid,grid},DifferenceOrder->#][DifferentiationMatrix] &;mat = dm@Pseudospectral;mattest=dm[points-d];Union@Flatten@Chop@Normal[mat-mattest]. (And whenDifferenceOrderis an option forNDSolve`, this should be counted. ) – xzczd Jun 24 '20 at 01:41

score 4 · Answer 2 · edited Sep 08 '19 at 06:47

4

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}}}]]

edited Sep 08 '19 at 06:47

xzczd

65,995
9
163
468

answered Sep 06 '19 at 11:30

user21

39,710
8
110
167

Side note, there's a hidden syntax for setting Coordinates: {x, mygrid} // Flatten. It's accidentally found here: https://mathematica.stackexchange.com/a/240406/1871 – xzczd Feb 22 '21 at 04:36
@xzczd, look at that. I did not know that. This was implemented before I joined the company... Good find! – user21 Feb 22 '21 at 10:39

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

2 Answers2

Linked