How to find the differentiation matrix based on FEM?

Question

Assume that I have a non-equidistant grid of $n$ nodes, as follows:

ClearAll["Global`*"];
n = 10;
SeedRandom[123];
nx = Sort@RandomReal[{-1, 6}, n]

If I want to compute the differentiation matrix in Mathematica, I simply use the command NDSolve`FiniteDifferenceDerivative as follows:

opt = "DifferenceOrder" -> 2;
dudx = NDSolve`FiniteDifferenceDerivative[1, nx, opt]["DifferentiationMatrix"];
dudx // MatrixForm

This computes the weights of an FD method which is very useful when solving a differential equation.

Is there a way to compute a similar matrix based on the finite element method (FEM), I mean with FEM accuracy? Maybe, a computation of a stiffness matrix could help?

Can we test and compare the performance of such two matrices (based on FD and FEM) for approximating a derivative of a function or a simple differential equation?

I would be thankful if someone write some hints to derive such a matrix with FEM methodology.

Answer 1

Here is a code I wrote a while back to get you started:

Needs["NDSolve`FEM`"]
FiniteElementDerivative[order : {__Integer}, mesh_ElementMesh] /; 
  1 <= Length[order] <= 3 := 
 Block[{dim, nr, vd, sd, mdata, ccoef, pos, dcoef, cdata},
  dim = Length[order];
  nr = ToNumericalRegion[mesh];
  vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, 
      Table[Unique[X], {dim}]}];
  sd = NDSolve`SolutionData[{"Space"} -> {nr}];
  mdata = InitializePDEMethodData[vd, sd];

  ccoef = ConstantArray[0, dim];
  pos = Flatten[Position[order, 1]];
  ccoef[[pos]] = 1;

  dcoef = ConstantArray[0, dim];
  pos = Flatten[Position[order, 2]];
  dcoef[[pos]] = 1;
  dcoef = DiagonalMatrix[dcoef];

  (* "Pure ConvectionCoefficients" will trigger a warning *)
  Quiet[
   cdata = InitializePDECoefficients[vd, sd
     , "DiffusionCoefficients" -> {{dcoef}}
     , "ConvectionCoefficients" -> {{ccoef}}
     ], {InitializePDECoefficients::femcscd}];

  DiscretizePDE[cdata, mdata, sd]
  ]

Examples:

mesh = ToElementMesh[Line[{{0}, {20}}], "MeshOrder" -> 1];
{dXmatFEM, d2XmatFEM} = 
  FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{1}, {2}};

d2XmatFEM // Normal

Compare with a FDM matrix:

ng = Flatten[mesh["Coordinates"]];
{dXmatFDM, d2XmatFDM} = 
  NDSolve`FiniteDifferenceDerivative[#, {ng}, "DifferenceOrder" -> 1][
     "DifferentiationMatrix"] & /@ {{1}, {2}};
d2XmatFDM // Normal

Note that the matrices scale differently (you will see this best when you change the length of the region) and that boundaries are treated differently.

In 2D:

mesh = ToElementMesh[Rectangle[{0, 0}, {1, 1}]];
{dXmat, dYmat, d2Xmat, d2Ymat, d2XYmat, lap} = 
  FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{1, 
     0}, {0, 1}, {2, 0}, {0, 2}, {1, 1}, {2, 2}};

d2Xmat // MatrixPlot

The operators are consistent with themselves:

SparseArray[d2Xmat + d2Ymat - lap] // Norm
2.14093742815996`*^-15