# Poor convergence for first order PDE using FEM (Laplace's equation in a disk)

Consider Laplace's eqn in the unit disk, r^2 D[u[r, t], {r, 2}] + r D[u[r, t], {r, 1}] + D[u[r, t], {t, 2}] == 0. The general analytic solution is easy to work out, and we can target the n=1 solution (i.e. u = r Cos[t]) using Dirichlet boundary conditions. The problem is easily solved with the finite element method:

uf1 = NDSolveValue[{r^2 D[u[r, t], {r, 2}] + r D[u[r, t], {r, 1}] +
D[u[r, t], {t, 2}] == 0,
DirichletCondition[u[r, t] == Cos[t], r == 1 && t != 2 Pi],
PeriodicBoundaryCondition[u[r, t], t == 2 \[Pi],
Function[x, x - {0, 2 \[Pi]}]]}, u, {r, 0, 1}, {t, 0, 2 Pi}];


Now let's try an alternative way to solve the same problem, using a first order PDE, D[u[r, t], r] == -I/r D[u[r, t], t], for a complex scalar field u. By separating u into its real and imaginary components, one can show that each satisfy the above Laplace's equation. We attempt to solve the PDE numerically, matching the Dirichlet condition above (for the real part):

uf2 = NDSolveValue[{D[u[r, t], r] == -I/r D[u[r, t], t],
DirichletCondition[u[r, t] == Exp[I t], r == 1 && t != 2 Pi],
PeriodicBoundaryCondition[u[r, t], t == 2 \[Pi],
Function[x, x - {0, 2 \[Pi]}]]}, u, {r, 0, 1}, {t, 0, 2 Pi}];


First things look OK: both solutions seem to approximate the exact target solution:

uA[r_, t_] = r Cos[t];
{Plot3D[Re[uf1[r, t]] - uA[r, t], {r, 0, 1}, {t, 0, 2 Pi}], Plot3D[Re[uf2[r, t]] - uA[r, t], {r, 0, 1}, {t, 0, 2 Pi}]}


but the conventional method is much better.

The question is: why does the 1st order method perform so poorly compared with the conventional (2nd order) method?

For the first order problem, Mathematica issues the usual warning about the equation being convection dominated, but can't we ignore this because we know the solutions are nice and smooth? (Btw, artificial diffusion doesn't help.) Is there an issue with applying FEM to first order PDEs, even if they are known to be elliptic?

A scan over MaxCellMeasure shows how serious the problem is:

errList1 = {}; errList2 = {};
For[n = 3, n <= 10, n++,
mcm = 1./2^n;
mesh = ToElementMesh[FullRegion[2], {{0, 1}, {0, 2 Pi}},
MaxCellMeasure -> mcm];
Quiet[
uf1 = NDSolveValue[{r^2 D[u[r, t], {r, 2}] + r D[u[r, t], {r, 1}] +
D[u[r, t], {t, 2}] == 0,
DirichletCondition[u[r, t] == Cos[t], r == 1 && t != 2 Pi],
PeriodicBoundaryCondition[u[r, t], t == 2 \[Pi],
Function[x, x - {0, 2 \[Pi]}]]}, u, {r, t} \[Element] mesh];
uf2 = NDSolveValue[{D[u[r, t], r] == -I/r D[u[r, t], t],
DirichletCondition[u[r, t] == Exp[I t], r == 1 && t != 2 Pi],
PeriodicBoundaryCondition[u[r, t], t == 2 \[Pi],
Function[x, x - {0, 2 \[Pi]}]]}, u, {r, t} \[Element] mesh];
];
err1 = Sqrt[
Mean[Table[{ri, ti} = c; (Re[uf1[ri, ti]] - uA[ri, ti])^2, {c,
mesh["Coordinates"]}]]];
err2 = Sqrt[
Mean[Table[{ri, ti} = c; (Re[uf2[ri, ti]] - uA[ri, ti])^2, {c,
mesh["Coordinates"]}]]];
AppendTo[errList1, {mcm, err1}];
AppendTo[errList2, {mcm, err2}];
];
ListLogLogPlot[{errList2, errList1}]


Although the error of the first solution goes clearly as mcm^2, the error of the second solution does not even convincingly go to zero with mcm.

A final motivating comment: Although the above problem is easily evaded (though not understood) by using the usual second order formulation, one may be faced with a first order PDE that cannot be easily transformed into the standard 2nd order form. It would be useful to know how reliably solve an elliptic first order system.

I can offer a partial answer to this. There are two issues I found.

(1) The polar coordinates were problematic at high resolution (I suppose the coordinate singularity at the origin might have been causing a problem), which did not allow me to go to very small mesh measure. By changing to cartesian coordinates, I was able to push to very small resolution and get more convincing convergence with the first order PDE:

m = 3; uA[x_, y_] = (x + I y)^m;
errList = {};

For[n = 5, n <= 14, n++,
mcm = 1./2^n;
meshxy = ToElementMesh[Disk[], MaxCellMeasure -> mcm];
Quiet[
uf = NDSolveValue[{(x - I y) D[u[x, y], x] + (y + I x) D[u[x, y],
y] == 0,
DirichletCondition[u[x, y] == uA[x, y], x^2 + y^2 == 1]}, u,
Element[{x, y}, meshxy]];
];
err = Sqrt[
Mean[Table[{xi, yi} = c; (Norm[uf[xi, yi] - uA[xi, yi]])^2, {c,
meshxy["Coordinates"]}]]];
AppendTo[errList, {mcm, err}];
];
ListLogLogPlot[errList]


The scaling of the error is only first order in the cell measure, but it is more convincing now.

(2) Poking around in FEM references, I noticed that first order PDEs have a reputation for being poorly handled by "standard" finite element methods, and some texts recommended a "least squares" approach (e.g. Ern, Guermond. Theory and Practice of Finite Elements (2004) chapter 5). This is equivalent to solving a related higher order PDE, but the format follows rather intuitively from the original PDE:

b = {x - I y, y + I x}; bs = {x + I y, y - I x};
PDE = -Div[bs b.Grad[u[x, y], {x, y}], {x, y}] == 0;


Now we do the error scan as before:

errList = {};
For[m = 5, m <= 14, m++, mcm = 1./2^m;
meshxy = ToElementMesh[Disk[], MaxCellMeasure -> mcm];
uf = NDSolveValue[{PDE,
DirichletCondition[u[x, y] == uA[x, y], x^2 + y^2 == 1]}, u,
Element[{x, y}, meshxy]];
err = Sqrt[
Mean[Table[{xi, yi} = c; (Norm[uf[xi, yi] - uA[xi, yi]])^2, {c,
meshxy["Coordinates"]}]]];
AppendTo[errList, {mcm, err}];
];
ListLogLogPlot[errList]


This error seems to scale as mcm^1.5.

EDIT

It may be possible to replicate the above performance, using FEM programming and post-processing, instead of solving the second order equation. Below is one idea that gives a slight improvement over the naive solution of the first order PDE (but does not do as well as solving the second order equation):

GetSystemMatricesLS[PDE_, iVars_, dVars_, mesh_] := Module[{nr, sd, vd, state, femdata, initBCs, initCoeffs, methodData, discretePDE, stiffness, load, LSstiffness, LSload, discreteBCs},
nr = ToNumericalRegion[mesh];
sd = NDSolveSolutionData[{"Space" -> nr}];
vd = NDSolveVariableData[{"DependentVariables" -> dVars,
"Space" -> iVars}];
methodData = InitializePDEMethodData[vd, sd];
Quiet[{state} =
NDSolveProcessEquations[PDE, dVars, iVars \[Element] mesh,
Method -> {"FiniteElement"}]];
femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
initCoeffs = femdata["PDECoefficientData"];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
stiffness = discretePDE["StiffnessMatrix"];
discreteBCs =
DiscretizeBoundaryConditions[initBCs, methodData, sd];
LSstiffness = ConjugateTranspose[stiffness].stiffness;
]


Note how the transformation to a symmetrized stiffness matrix is performed before deploying the boundary conditions.

k = 3; uA[x_, y_] = (x + I y)^k;

errList = {};
For[n = 5, n <= 14, n++,
mcm = 1./2^n;
PrintTemporary["n: ", n];
meshxy = ToElementMesh[Disk[], MaxCellMeasure -> mcm];
GetSystemMatricesLS[{(x - I y) D[u[x, y], x] + (y + I x) D[u[x, y],
y] == 0, DirichletCondition[u[x, y] == uA[x, y], True]}, {x,
y}, {u}, meshxy];

• (+1) This is very interesting. I have fixed a syntax error in your second part of the code and removed the Quiet. If you could post the FEM programming method that you have I have an idea of how would could get this to work from NDSolve. We could specify a custom LinearSolve for this. I tried this but you not get it to work with the LinearSolve command you provided. If I could have a look at your code I might be able to figure this out how to do this. Thanks. – user21 Nov 28 '19 at 8:09
• Maybe I miss read your answer. I was under the impression that the second error plot could be achieved with LinearSolve[ConjugateTranspose[stiffness].stiffness, ConjugateTranspose[stiffness].load], is that correct? I put the code in an "answer" that will delete once you read it. – user21 Nov 28 '19 at 20:37
• Hello again. You are right, I did expect the error plots to be similar, but they clearly are not. Appears that I just jumped to the wrong conclusion, based on an encouraging plots, and should been more careful. After experimenting a bit, I did manage to find some improvement with a post-processing technique, by delaying the application of DeployBoundaryConditions` until after I have applied the transformation into a Hermitian stiffness matrix. I will update the "Solution" with that code. – Will.Mo Nov 29 '19 at 12:48