# NDSolveFiniteDifferenceDerivative gives wrong result when the precision is not MachinePrecision

Bug introduced in 8 or earlier and persisting through 11.0.1, fixed in 11.3 or earlier.

I want to get a pseudospectral differentiation matrix by NDSolveFiniteDifferenceDerivative.

nx = 8;  (* Grid size *)
x = N[Table[-Cos[(j \[Pi])/nx], {j, 0, nx}]];
(* -> {-1., -0.92388, -0.707107, -0.382683, 0., 0.382683, 0.707107, 0.92388, 1.} *)


The differentiation matrix is

dx = NDSolveFiniteDifferenceDerivative[1, x, DifferenceOrder -> "Pseudospectral"]["DifferentiationMatrix"]


We can check that dx.x is 1.

dx.x
(* -> {1., 1., 1., 1., 1., 1., 1., 1., 1.} *)


However, if we change the precision to 20 (arbitrary precision), the result is obviously wrong.

nx = 8;
x1 = N[Table[-Cos[(j \[Pi])/nx], {j, 0, nx}], 20]
dx1 = NDSolveFiniteDifferenceDerivative[1, x1, DifferenceOrder -> "Pseudospectral"]["DifferentiationMatrix"]


dx1.x1 is supposed to be 1, but it is far from 1.

dx1.x1
(* -> {22.50000000000000000, 2.649846074667681225, 0.535533905932737622, 1.323575230741387329, 0.500000000000000000, 2.451028024842430524, -6.535533905932737622, 125.57555066974850092, 252.50000000000000000} *)


Some elements of dx1-dx are far from zero (though most of its elements are close to zero).

Similar things happen for finite difference and periodic differentiation matrices. Is there a way to increase the precision and get the correct result?

• The last column seems in error to me. Yes? – Michael E2 Jan 16 '16 at 20:28
• Yes, only the last column of dx1-dx is far from zero. – renphysics Jan 16 '16 at 20:33
• I filed this as a bug. Thanks for reporting. – user21 Jan 17 '16 at 16:46
• Oh, actually I always found that "DifferenceOrder" -> "Pseudospectral" seems not to work properly when a non-MachinePrecision WorkingPrecision is set, maybe this is the true reason? – xzczd Oct 2 '16 at 15:28

The last column seems in error. Here's a workaround for the sample problem, although it does not fix NDSolveFiniteDifferenceDerivative:

dx2 = dx1;
dx2[[All, -1]] = -Reverse@dx1[[All, 1]]
(*
{-0.50000000000000000000, 0.25989153247414500869,
-0.29289321881345247560, 0.36161567304292239214,
-0.50000000000000000000, 0.80995720221088751026,
-1.7071067811865475244, 6.5685355922720450889, 21.500000000000000000}
*)

dx2.x1
(*
{1.00000000000000000, 1.000000000000000000, 1.000000000000000000,
1.000000000000000000, 1.000000000000000000, 1.000000000000000000,
1.000000000000000000, 1.000000000000000000, 1.00000000000000000}
*)


Hope that helps.

Update 1: Interestingly, this works:

NDSolveFiniteDifferenceDerivative[1, x1, x1, DifferenceOrder -> "Pseudospectral"]
(*  {1.00000000000000000,..., 1.00000000000000000}  *)


So does this:

dx1FN = NDSolveFiniteDifferenceDerivative[1, x1, DifferenceOrder -> "Pseudospectral"]
dx1FN[x1]


Update 2: Here's another workaround to get the differentiation matrix:

Transpose[dx1FN /@ IdentityMatrix[Length@x1]]


Update 3: The "DifferentiationMatrix" is computed internally basically by calling

NDSolveFiniteDifferenceDerivative[Derivative[1, 0], {x1, x1},
DifferenceOrder -> "Pseudospectral"][IdentityMatrix[Length@x1]]


This produces the erroneous matrix and should probably be considered a bug. Strangely, it computes the $(0,1)$ derivative correctly. So here is yet another workaround:

Transpose[
NDSolveFiniteDifferenceDerivative[Derivative[0, 1], {x1, x1},
DifferenceOrder -> "Pseudospectral"][IdentityMatrix[Length@x1]]]
`