13
$\begingroup$

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 NDSolve`FiniteDifferenceDerivative.

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 = NDSolve`FiniteDifferenceDerivative[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 = NDSolve`FiniteDifferenceDerivative[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?

$\endgroup$
  • $\begingroup$ The last column seems in error to me. Yes? $\endgroup$ – Michael E2 Jan 16 '16 at 20:28
  • $\begingroup$ Yes, only the last column of dx1-dx is far from zero. $\endgroup$ – renphysics Jan 16 '16 at 20:33
  • 3
    $\begingroup$ I filed this as a bug. Thanks for reporting. $\endgroup$ – user21 Jan 17 '16 at 16:46
  • $\begingroup$ 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? $\endgroup$ – xzczd Oct 2 '16 at 15:28
13
$\begingroup$

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

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:

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

So does this:

dx1FN = NDSolve`FiniteDifferenceDerivative[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

NDSolve`FiniteDifferenceDerivative[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[
 NDSolve`FiniteDifferenceDerivative[Derivative[0, 1], {x1, x1}, 
   DifferenceOrder -> "Pseudospectral"][IdentityMatrix[Length@x1]]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.