9
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ This FEM tutorial shows how to get the system matrices. $\endgroup$ – Michael E2 May 21 '18 at 21:54
  • $\begingroup$ To compare the results of the two different ways, it would be good if we test them on a differential equation solving. $nx$ could be changed to a better distribution of the nodes (not random) and then compare practically. $\endgroup$ – M.J.2 May 21 '18 at 22:11
8
$\begingroup$

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

enter image description here

The operators are consistent with themselves:

SparseArray[d2Xmat + d2Ymat - lap] // Norm
2.14093742815996`*^-15
$\endgroup$
3
  • $\begingroup$ A very good response. However, to complete the discussion I have several questions. 1. Why the incorporation of the boundaries are different for these two derivative matrices (FDM and FEM)? For the derivative matrix of the FDM, we basically change the first and last rows by the boundaries (when solving a time-dependent 1D PDE). 2. Can we use the KroneckerProduct[] to obtain the derivative matrices for the 2D and 3D cases when using FEM? Since, for the derivative matrix of the FDM, it is possible to find the derivative matrix in higher dimensions using this? $\endgroup$ – Fazlollah May 22 '18 at 15:36
  • $\begingroup$ In addition, 3. Does the derivative matrix of the FEM, rely on the structure of the mesh (in the 1D case)? Because, I think whatever the structure is (equidistant or non-equidistant), the derivative matrix would be the same (for the first derivative-approximation)? $\endgroup$ – Fazlollah May 22 '18 at 15:38
  • $\begingroup$ @Fazlollah, 1) I think the theory is a bit different for FEM/FDM. 2) How the boundaries are changed in the system matrices depends on the boundaries you want to model. For for FEM not changing anything means Neumann 0 BCs. To the best of my knowledge KroneckerProduct can not be used for higher order FEM derivation. 3) You can try that out: Make an ElementMesh with non equidistant coords. $\endgroup$ – user21 May 23 '18 at 5:36

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.