# FEM doesn't output exactly the same solution for exactly the same code?

Consider this toy example:

sol :=
NDSolveValue[{Laplacian[u[x, y], {x, y}] + 1 == 0,
DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}]

sol == sol
(* True *)

sol - sol
(* - InterpolatingFunction[…] + InterpolatingFunction[…] *)

sol === sol
(* False *)

Table[sol[0.6, 0.6] - sol[0.6, 0.6], {10}]
(* {3.05311*10^-16, -1.249*10^-16, -9.71445*10^-17,
-4.85723*10^-16, 9.71445*10^-17, 2.77556*10^-17,
-5.55112*10^-17, 0., 1.80411*10^-16, 2.498*10^-16}  *)


Tested in v12.3 and v13.2.

What's happening here?

Is it possible to force NDSolveValue to produce exactly the same output for every execution?

• It is not so clear what do you try to test? Nor any numerical parameters are fixed. Can you compare sol on the same mesh? Commented Jan 6, 2023 at 6:10
• @alex I'm testing the equivalence of different PDE systems. (To be more specific, I'm exploring this question: mathematica.stackexchange.com/q/278215/1871 ) I've tried generating mesh beforehand, the error is still there. Commented Jan 6, 2023 at 6:20
• I think the problem comes from "SameQ". It is very specific what it considers to be True. Commented Jan 6, 2023 at 8:40

You are right to expected that a deterministic algorithm should produce the exact same result every time it is run. And solving the linear equations of this elliptic PDE should be deterministic. Or not?

In fact, there are several algorithms that use some randomness. Mesh generation is such a thing as user21 mentioned.

I'd like to point to another source of randomness: the linear solver that finally solves the linear equations. Or more precisely the algorithm that permutes rows and columns of the sparse system matrix. Here the aim is to find a permutation such that Gaussian elimination introduces as little fill-in as possible (and some other subtle side conditions, e.g. that the induced elimination tree is well balanced and not too deep ). This is a NP hard problem. So one has to settle with heuristics. I don't know which exact library is employed here. But I would not be surprised if it employs some randomization. That would imply that from run to run of the solver the Gauss algorithm eliminates the rows of the system matrix in different orders -- which would not matter in exact arithmetic but introduces some rounding errors in floating point arithmetic.

• @xzczd what your all think, would you expect this type of information in the documentation?, If yes, where would you place it? Maybe a possible issue for NDSolve? Commented Jan 6, 2023 at 9:06
• This seems to be the correct direction. See my answer. Commented Jan 6, 2023 at 9:19
• @user21 I think possible issue of NDSolve is the place to go. Commented Jan 6, 2023 at 9:22
• @xzczd, I have added something to the LinearSolve and NDSolve possible issues sections. Commented Jan 6, 2023 at 13:31
• @Henrik How can you explain my example? :) Commented Jan 6, 2023 at 15:32

Inspired by Henrik's answer, I think I find the root of problem: "Pardiso" method of LinearSolve is the culprit. According to the tutorial NDSolve Options for Finite Elements:

The finite element method uses the efficient direct PARDISO solver as the default linear solver.

Let's test Pardiso:

n = 1000; m =
SparseArray[{{i_, i_} -> -2., {i_, j_} /; Abs[i - j] == 1 -> 1.}, {n, n}];
b = Table[0., {n}]; b[[99]] = 1.;

sollinear := LinearSolve[m, b, Method -> "Pardiso"]

sollinear - sollinear // Abs // Max
(* 2.13163*10^-12 *)


It's worth noting that n should be large enough. For small system e.g. n = 100, the randomness becomes 0.. This explains why Alex's method works. (MaxCellMeasure -> 10^-2 actually results in a mesh coarser than default. )

If we turn to following methods for LinearSolve, deterministic solution will be produced (in reasonable time):

solmethod[method_] :=
NDSolveValue[{Laplacian[u[x, y], {x, y}] + 1 == 0,
DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"PDESolveOptions" -> {"LinearSolver" -> {Automatic, Method -> method}}}}]

Table[
solmethod[#] === solmethod[#] & /@ {"Multifrontal", "Krylov",
"IterativeRefinement", "Banded"}, {10}] // Flatten // Union
(* {True} *)


Still, this isn't the end: what if the Pardiso method is necessary for certain equation solving? Is it possible to adjust option of Pardiso method to remove its randomness?

• Aha. Yes. I recall that Pardiso had some option to turn on CNR mode (CNR = conditional numerical reproducibility). @user21 One can switch it on by modifying some of the infamous iparm parameters when initializing the Pardiso solver. Have a look at this page of Pardiso's documentation and watch out for "CNR". Commented Jan 6, 2023 at 9:37
• Alas, that will switch off parallelization of computing the nested dissection ordering. And doubt that anybody really wants to do that... Commented Jan 6, 2023 at 9:38
• @HenrikSchumacher, thanks, I had a look. However exposing that option is not going to happen. The amount of work and unexpected trouble this will produce is not warranted for such a 'cosmetic' improvement IMHO, Commented Jan 6, 2023 at 13:29
• "For small system e.g. n = 100, the randomness becomes 0." That makes sense. The Pardiso doc page states that randomness is introduced by only for parallel computation of the matrix reordering. And parallelization is probably truended off for too small matrices. Commented Jan 7, 2023 at 6:52
• By the way, I think that Pardiso uses Metis for computing the reordering. I skimmed a bit through the documentation of the multithreaded version MT-Metis. Could be that the "randomness" is entirely due to race conditions. (I was somewhat chocked to read that!) The Metis authors state that this were not an issue because they also make sure that the reordering is valid. Nonetheless, the race conditions make the library quite nondeterministic. Commented Jan 7, 2023 at 6:55

Contra example

Needs["NDSolveFEM"]
mesh = ToElementMesh[FullRegion[2], {{0, 1}, {0, 1}},
MaxCellMeasure -> 10^-2]

(*Out[]= ElementMesh[{{0., 1.}, {0., 1.}}, {QuadElement["<" 100 ">"]}]*)

Do[
sol[i] =
NDSolveValue[{Laplacian[u[x, y], {x, y}] + 1 == 0,
DirichletCondition[u[x, y] == 0, True]}, u,
Element[{x, y}, mesh]];, {i, 10}]

Table[sol[1][6/10, 6/10] - sol[i][6/10, 6/10], {i, 10}]

{0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

• If MaxCellMeasure -> 10^-2 is removed, the random error is there again. Might be because the error is smaller than $MachineEpsilon for very fine mesh. Still, this is an answer to my last question. (+1) Commented Jan 6, 2023 at 6:48 • I am not a 100% sure but I believe that the mesh generation has some degree of randomness to seed the algorithms. Commented Jan 6, 2023 at 6:58 • "Might be because the error is smaller than$MachineEpsilon for very fine mesh." I retract this line, because with MaxCellMeasure -> .001 the randomness shows up again, so MaxCellMeasure -> 10^-2 is actually too coarse… Commented Jan 7, 2023 at 2:39
• @xzczd Nevertheless we have same answer for every run. Commented Jan 9, 2023 at 12:45
• You mean for MaxCellMeasure -> .001 the error is constant? It's not the case on my laptop (Win 10) and Wolfram Cloud (Linux): wolframcloud.com/obj/a4938367-256c-4772-a0d2-11baa90c095a Commented Jan 9, 2023 at 14:01