# Inverting differential equation using finite element methods

tl;dr; How to use FEM tools to extract models needed to invert PDEs.

Context

In astrophysics, one is interested in so-called 'cosmic archeology' which involves recovering the origin of a given observation, while modelling its evolution. The idea is to be able to understand what may have caused in the past an given signature found in the data. For instance, can we explain the observed vertical velocity distribution of stars above and below the galactic disc seen by the Gaia spacecraft in terms of past satellites having hit our Milky way. Example

As a test example let us consider a 1D diffusion equation sourced by a finite set of heat sources

 source[x_, t_] =
3 BSplineBasis[3, t 4] BSplineBasis[3, (x - 1/8) 4] +
2 BSplineBasis[3, -2 + t 4] BSplineBasis[3, (x - 5/8) 4] +
BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 2/8) 4] +
BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 1/8) 4] +
BSplineBasis[3, -1/2 + t 4] BSplineBasis[3, (x - 4/8) 4] +
3/2 BSplineBasis[3, -3 + t 4] BSplineBasis[3, (x - 1/8) 4];

ContourPlot[source[x, t], {x, 0, 1}, {t, 0, 1}, PlotPoints -> 36,
Exclusions -> None, PlotRange -> All,
AspectRatio -> 1, Contours -> 10] The diffusion diagram will obey

    sol0 = NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] ==
source[x, t],
f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0}, f, {x, 0, 1}, {t, 0, 2}];

ContourPlot[sol0[x, t], {x, 0, 1}, {t, 0, 1}, FrameLabel -> {x, t},
AspectRatio -> 1, PlotRange -> All, Contours -> 30, PlotPoints -> 50] Here I have assumed arbitrarily that the edges of [0,1] did not let heat diffuse. I also assumed that initially there was no heat.

Let me first generate the corresponding data set of positing and time $$(x,t)$$ for later use

data = N[Flatten[
Table[{x, t, sol0[x, t]}, {x, 0, 1, 1/32}, {t, 0, 1, 1/32}], 1]];


My purpose is to invert this data set to recover the source of heating.

I other words, can I recover the first plot from the second one, if I assume I know how the heat source diffuse?

Attempt

I can define a set of spline functions which cover the $$(x,t)$$ space as follow:

nn = 16;
knots = Flatten[{{0, 0}, (Range[0, nn]/nn), {1, 1}}];

basis0 = Flatten@
Table[BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j,
t], {i, 0, nn}, {j, 0, nn}];


For instance, the 60th function obeys:

Plot3D[basis0[], {x, 0, 1}, {t, 0, 1}, PlotRange -> All] The image of this basis satisfies

basis = Flatten@
Table[NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] ==
BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, t],
f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0},
f[x, t], {x, 0, 1}, {t, 0, 1}], {i, 0, nn}, {j, 0, nn}];

Plot3D[basis[], {x, 0, 1}, {t, 0, 1}, PlotRange -> All] From this basis and the data I can generate the linear model a0 and a relating resp. the spline coefficients to the source map and the spline coefficients heat map:

ff = Function[{x, t}, basis0 // Evaluate];
a0 = ff @@ # & /@ (Most /@ data0);


and

ff = Function[{x, t}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data); Let me first check that I can fit the source map with the splines:

fit0[x_, t_] =
basis0.LinearSolve[Transpose[a0].a0, Transpose[a0].(Last /@ data0)];
ContourPlot[fit0[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20,
PlotRange -> All] Similarly, I can define an (isotropic) penalty corresponding to $$\int |\Delta T|^2 dx dt$$ as

ff = Function[{x, t}, D[basis0, x, x] + D[basis0, t, t] // Evaluate];
s0 = ff @@ # & /@ (Most /@ data0);
pen = Transpose[s0].s0; pen /= Max[Flatten[Abs[pen]]]; A solution to the inverse problem then follows simply from inverting a with a small roughness penalty as

sol[x_, t_] =
basis0.LinearSolve[Transpose[a].a + 10^-7 pen,
Transpose[a].(Last /@ data)];

ContourPlot[sol[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20,
PlotRange -> All] Question

I am fairly certain my present implementation is effectively redundant with the way NDSolve can actually solves the differential equation using Finite Element methods. Hence my question:

Could one make use of the actual solver in NDSolve to formulate the inverse problem?

In other words, can we extract from the FEM toolkit FEM_a,FEM_source and FEM_solution and FEM_grid so that

   FEM_solution =  FEM_a  FEM_source


where 'FEM_' stands for as sampled by the underlying mesh,FEM_grid of the FEM toolkit.

This would be of interest in terms of efficiency, but also in order to address more complex and realistic inverse problems?

For Gaia data, the diffusion is in fact occurring in 3D and is anisotropic, so a robust and efficient formulation would help!

Technically I believe FEM have access to both a and pen so it would be great to be able to access them for the sake of solving the inverse problem.

I am guessing that this link would be a good starting point?

Comment

Note that the above implementation is partially incorrect at the top edge, because most spline basis elements are required to be zero on the boundary, whereas the correct solution should have an outgoing flux condition. This is something the FEM would handle naturally when the boundary conditions are taken care of.

  Plot[{sol[x, 1], sol0[x, 1]}, {x, 0, 1}] Note that in astronomy we unfortunately don't have access to the full diffusion diagram but only typically to a given snapshot (i.e. data on a line at fixed time, and/or possibly the time derivative on that line). So we can only extrapolate in the past up to some quite limited time horizon.

Complement 1: 1+1D code

source[x_, t_] =
3 BSplineBasis[3, t 4] BSplineBasis[3, (x - 1/8) 4] +
2 BSplineBasis[3, -2 + t 4] BSplineBasis[3, (x - 5/8) 4] +
BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 2/8) 4] +
BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 1/8) 4] +
BSplineBasis[3, -1/2 + t 4] BSplineBasis[3, (x - 4/8) 4] +
3/2 BSplineBasis[3, -3 + t 4] BSplineBasis[3, (x - 1/8) 4];
sol0 = NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] ==
source[x, t],
f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0},
f, {x, 0, 1}, {t, 0, 2}];
nn = 16; knots = Flatten[{{0, 0}, (Range[0, nn]/nn), {1, 1}}];
basis0 = Flatten@
Table[BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j,
t], {i, 0, nn}, {j, 0, nn}];
basis = Flatten@
Table[NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] ==
BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, t],
f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0},
f[x, t], {x, 0, 1}, {t, 0, 1}], {i, 0, nn}, {j, 0, nn}];
data = N[Flatten[
Table[{x, t, sol0[x, t]}, {x, 0, 1, 1/32}, {t, 0, 1, 1/32}], 1]];
data0 = N[
Flatten[Table[{x, t, source[x, t]}, {x, 0, 1, 1/32}, {t, 0, 1,
1/32}], 1]];
ff = Function[{x, t}, basis0 // Evaluate];
a0 = ff @@ # & /@ (Most /@ data0);
ff = Function[{x, t}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data);
fit0[x_, t_] =
basis0.LinearSolve[Transpose[a0].a0,
Transpose[a0].(Last /@ data0)];
fit[x_, t_] =
basis.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)];
ff = Function[{x, t}, D[basis0, x, x] + D[basis0, t, t] // Evaluate];
s0 = ff @@ # & /@ (Most /@ data0);
pen = Transpose[s0].s0; pen /= Max[Flatten[Abs[pen]]];
sol[x_, t_] =
basis0.LinearSolve[Transpose[a].a + 10^-7 pen,
Transpose[a].(Last /@ data)];
ContourPlot[source[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20,
PlotRange -> All,Exclusions -> None]
ContourPlot[sol[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20,
PlotRange -> All]


Complement 2: 2+1D codes

For the sake of completeness and to demonstrate why a more efficient implementation is needed here is the code for 2D diffusion without FEM (which for n=16 would take a white to run!).

source[x_, y_, t_] = BSplineBasis[3, t ] BSplineBasis[3, x]*
BSplineBasis[3, y]
sol0 = NDSolveValue[{D[f[x, y, t], t] - 1/4 D[f[x, y, t], x, x] -
1/4 D[f[x, y, t], y, y] == source[x, y, t], f[x, y, 0] == 0,
DirichletCondition[f[x, y, t] == 0, True]}, f, {x, 0, 1}, {y, 0, 1}, {t, 0, 1}]
nn = 2;knots = Flatten[{{0, 0}, (Range[0, nn]/nn), {1, 1}}];
basis0 = Flatten@
Table[BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, y]
BSplineBasis[{3, knots}, k, t], {i, 0, nn}, {j, 0, nn}, {k, 0, nn}];
basis = Flatten@(Table[
ParallelTable[
NDSolveValue[{D[f[x, y, t], t] - 1/4 D[f[x, y, t], x, x] -
1/4 D[f[x, y, t], y, y] ==
BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j,
y] BSplineBasis[{3, knots}, k, t], f[x, y, 0] == 0,
DirichletCondition[f[x, y, t] == 0, True]},
f[x, y, t], {x, 0, 1}, {y, 0, 1}, {t, 0, 1}], {j, 0, nn}, {k,
0, nn}], {i, 0, nn}]);
data0 = N[Flatten[Table[{x, y, t, source[x, y, t]}, {x, 0, 1, 1/nn/2},
{y, 0,1, 1/nn/2}, {t, 0, 1, 1/nn/2}], 2]];
data = N[Flatten[
Table[{x, y, t, sol0[x, y, t]}, {x, 0, 1, 1/nn/2}, {y, 0, 1,
1/nn/2}, {t, 0, 1, 1/nn/2}], 2]];
ff = Function[{x, y, t}, basis // Evaluate];
a = ParallelMap[ff @@ # &, Most /@ data];
ff = Function[{x, y, t}, D[basis0, x, x] +
D[basis0, y, y] + D[basis0, t, t] // Evaluate];
s0 = ff @@ # & /@ (Most /@ data0);
pen = Transpose[s0].s0; pen /= Max[Flatten[Abs[pen]]];
sol[x_, y_, t_] =
basis0.LinearSolve[Transpose[a].a + 10^-9 pen,
Transpose[a].(Last /@ data)];
ContourPlot[sol[x, 1/2, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20,
PlotRange -> All]


Complement 3: Background

Let

$$\mathcal{L}\psi = \rho$$

represent a (linear) partial differential equation (possibly time dependant). I will assume that there exist a basis function over which I can project $$\psi$$, so that $$\psi(x)=\sum_n a_n \phi_n(x)\,,$$ where I also request that all $$\phi_n(x)$$ satisfy the boundary conditions of the partial differential equation ( $$x$$ can represent say $$\mathbf{r}$$ or $$(\mathbf{r},t)$$), i.e. the analysis is not necessary limited to stationary PDE). If I put this expansion into the PDE, multiply by $$\phi_p(x)$$ (or a Dirac function as a variant, see below) and integrate over $$x$$, I get formally $$\mathbf{L}\cdot \mathbf{\Phi} = \mathbf{P}\,,$$ where $$L_{ij}= \int d x \phi_i \mathcal{L} \phi_j$$, $$P_{i}= \int d x \phi_i \rho$$ and $${\Phi}_i= a_i$$.

I can then invert for $$\mathbf{\Phi}$$ as $$\mathbf{\Phi} =\mathbf{L}^{(-1)} \cdot\mathbf{P}\,,$$ where $$\mathbf{L}^{(-1)}$$ is the (possibly regularised) pseudo inverse of $$\mathbf L$$ (e.g. through least square). This is a possible method for solving PDEs. I am assuming (wrongly?) that linear FEM methods are a variant of this technique.

$$\psi = \mathcal{L}^{-1}\rho \,.$$ I can expand $$\rho$$ over a basis function,$$\rho=\sum_n a_n \rho_n$$ , project as previously and write eventually

$$\mathbf{P} =\mathbf{R}^{(-1)}\cdot \mathbf{\Phi}\,,$$ where $$\mathbf{R}^{(-1)}$$ is the (possibly regularised) pseudo inverse of $$\mathbf R$$, whose components are $$R_{ij}= \int d x \rho_i \mathcal{L}^{-1} \rho_j$$. In my code above I have implemented something closely related to the second method, relying on NDSolve (I use a Dirac function instead of $$\rho_i$$ to simply sample the measured $$\phi(x)$$). The reason is I am after recovering the source $$\mathbf P$$ given some knowledge of the response $$\mathbf \Phi$$.

My hope is that since FEM method solve the first problem they should have in store the tools to solve the second problem more efficiently?

• Oh finally here comes a serious question about inverse problem! I know little about the underlying theory, but am eager to learn something from the (possibly) upcoming answer(s)! Mar 22, 2020 at 3:55
• @xzczd If you are interested in inverse pbs I strongly recommend the web page of the late Albert Tarantola ipgp.jussieu.fr/%7Etarantola Mar 22, 2020 at 20:20
• Personally I learnt from that paper ipgp.jussieu.fr/%7Etarantola/Files/Professional/Papers_PDF/… Mar 22, 2020 at 20:34
• While I am interested I am not sure I understand the question. You want to write the above code with low level functions, in order to minimize computational time? I can certainly help in getting a low level FEM code of the ground but I am no expert in inverse problems. Some questions: What is the purpose of a0 it does not seem to be used. Also your "complement" code misses a definition of s0. Have you had a look at the FEM programming tutorial? Mar 23, 2020 at 6:18
• @chris, not entirely, I must admit. My feeling is that this is project that would require more thought and time on my side to fully understand. A commitment which I do not want to make right now. Sorry about that. Mar 25, 2020 at 11:41

I am going to show how to write this part of your post

source[t_, x_] =
3 BSplineBasis[3, t 4] BSplineBasis[3, (x - 1/8) 4] +
2 BSplineBasis[3, -2 + t 4] BSplineBasis[3, (x - 5/8) 4] +
BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 2/8) 4] +
BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 1/8) 4] +
BSplineBasis[3, -1/2 + t 4] BSplineBasis[3, (x - 4/8) 4] +
3/2 BSplineBasis[3, -3 + t 4] BSplineBasis[3, (x - 1/8) 4];

tEnd = 2;
AbsoluteTiming[
sol0 = NDSolveValue[{D[f[t, x], t] - 1/4 D[f[t, x], x, x] ==
source[t, x], f[0, x] == 0, f[t, 0] == 0, f[t, 1] == 0},
f, {x, 0, 1}, {t, 0, tEnd}
, Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}
];]
(* {0.337181, Null} *)


with the low level FEM functions. It's still not quite clear to me how you want to make use of this. More on this later. Note that I added a method option to tell NDSolve to actually make use of the FEM. Not all of the NDSolve calls you show actually make use of the FEM. But I think the method used is also not relevant.

To understand the code that follows I'd advise to read the FEMProgramming tutorial.

Set up the region, bcs, ics:

region = Line[{{0}, {1}}];
bcs = {DirichletCondition[u[t, x] == 0, True]};
initialConditionValue = 0.;
vd = NDSolveVariableData[{"DependentVariables" -> {u},
"Space" -> {x}, "Time" -> t}];

Needs["NDSolveFEM"]
nr = ToNumericalRegion[region];
sd = NDSolveSolutionData[{"Space" -> nr, "Time" -> 0.}];


Set up the PDE coefficients without the load term:

dim = RegionDimension[region];
initCoeffs =
InitializePDECoefficients[vd,
sd, {"DampingCoefficients" -> {{1}},
"DiffusionCoefficients" -> {{-1/4 IdentityMatrix[dim]}}}];


We omit the load term for now, as that is the term that is variable in your examples and we will take care of that later.

Initialize the BCs, method data and compute the stationary (time independent) discretization and boundary conditions of the PDE (without the load). These coefficients and discretizations are the same for all the PDEs you solve so we compute them only once.

initBCs = InitializeBoundaryConditions[vd, sd, {bcs}];
methodData = InitializePDEMethodData[vd, sd];

sdpde = DiscretizePDE[initCoeffs, methodData, sd, "Stationary"];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];


Now, we want to write a residual function for NDSolve to time integrate. At the same time we want the source to be variable.

makeResidualFunction[load_] := With[
InitializePDECoefficients[vd,
With[
discretizePDEResidual[t_?NumericQ, u_?VectorQ, dudt_?VectorQ] :=

NDSolveSetSolutionDataComponent[sd, "Time", t];
NDSolveSetSolutionDataComponent[sd, "DependentVariables", u];

{l, s, d, m} = sdpde["SystemMatrices"];

(* discretize and add the laod *)

DeployBoundaryConditions[{l, s, d}, sbcs];

d.dudt + s.u - l
]
]
]


This functions get the 'source' function and defines a residual function.

Generate the initial conditions with boundary conditions deployed.

init = Table[
initialConditionValue, {methodData["DegreesOfFreedom"]}];
init[[sbcs["DirichletRows"]]] = Flatten[sbcs["DirichletValues"]];


Get the sparsity pattern for the damping matrix for the time integration.

sparsity = sdpde["DampingMatrix"]["PatternArray"];


Set up the residual function.

makeResidualFunction[source[t, x]]


Do the time integration

AbsoluteTiming[
ufun = NDSolveValue[{
discretizePDEResidual[t, u[t], u'[ t]] == 0
, u == init}, u, {t, 0, tEnd}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {Automatic, Sparse -> sparsity}
(*,EvaluationMonitor\[RuleDelayed](monitor=Row[{"t = ",CForm[t]}])*)
, AccuracyGoal -> $$MachinePrecision/4, PrecisionGoal ->$$MachinePrecision/4
]
]

(* {0.429631,.... *)


As you see the time integration is somewhat slower from top level code.

Convert the result to an interpolating function:

ffun = ElementMeshInterpolation[{ufun["Coordinates"][],
methodData["ElementMesh"]}, Partition[ufun["ValuesOnGrid"], 1]]


Check that this is on the same order as the NDSolve result.

Plot3D[sol0[t, x] - ffun[t, x], {t, 0, tEnd}, {x, 0, 1},
PlotRange -> All]


Discussion:

I think you make an implicit assumption that is not correct. You assume that the matrix assembly process is the expensive part. But it's not. It's the actual time integration that you need to do many many times that is expensive. Precomputing the system matrices can probably save a little when the parallel computation is run but you can not make the time integration go away altogether.

• thanks very much. Yes I did wonder about the time integration per basis element. Mar 23, 2020 at 10:33

Let me try and begin to answer my own question, as it might inspire better answers. Here I will solve the Poisson equation as a test case using 0-splines.

Needs["NDSolveFEM"];
reg0 = Rectangle[{0, 0}, {1, 1}];
mesh0 = ToElementMesh[reg0, MaxCellMeasure -> 0.025, AccuracyGoal -> 1]


I can then extract the mesh elements

idx = mesh0["MeshElements"][[1, 1]];mesh0["Wireframe"] In order to define the density on each cell I need to find the convex hull of each cell

pol = Table[mesh0["Coordinates"][[ idx[[i]]]] // ConvexHullMesh, {i,Length[idx]}];


I can then use the function RegionMember to define the Indicator of that cell (as shown in this answer)

basis = Table[f[x_, y_] = Boole[ RegionMember[pol[[i]], {x, y}]];
NDSolveValue[{-Laplacian[u[x, y], {x, y}] == f[x, y]
+ NeumannValue[0, True] // Evaluate,DirichletCondition[u[x, y] == 0, True]},
u[x, y], {x, y} \[Element] mesh0],{i, Length[idx]}];


Then I can plot the basis

Plot3D[basis[[;; ;; 5]], {x, y} \[Element] mesh0,
PlotStyle -> Opacity[0.4], PlotRange -> All, PlotTheme -> "Mesh"] Of course the main point of using the mesh of the FEM is that it can be non trivial. For instance

Needs["NDSolveFEM"];
mesh0 = ToElementMesh[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]],
MaxCellMeasure -> 0.25, AccuracyGoal -> 1]; mesh0["Wireframe"] while the same code will work exactly unchanged

pol = Table[mesh0["Coordinates"][[ idx[[i]]]] // ConvexHullMesh, {i,Length[idx]}];
basis = Table[f[x_, y_] = Boole[ RegionMember[pol[[i]], {x, y}]];
NDSolveValue[{-Laplacian[u[x, y], {x, y}] == f[x, y] +
NeumannValue[0, True] // Evaluate,
DirichletCondition[u[x, y] == 0, True]},
u[x, y], {x, y} \[Element] mesh0],{i, Length[idx]}];


And once again

Plot3D[basis[[;; ;; 5]], {x, y} \[Element] mesh0,
PlotStyle -> Opacity[0.4], PlotRange -> All, PlotTheme -> "Mesh"]


Now the inverse problem is quite simple I find the FEM toolkit extremely useful in this context because building a basis function for non trivial geometry is ... non trivial, while the FEM package takes care of it all here.

This solution does not fully address my original question because the basis are 0-splines. Ideally cubic spline would be good too.

Proof of concept for inversion

Let's see how the basis can be used to fit some model. Let us start with a basis on the mesh

basis0 = Table[Boole[ RegionMember[pol[[i]], {x, y}]], {i,Length[idx]}];


source[x_, y_] = basis0[[{150, 170, 125}]].{2, 4, 5};
ContourPlot[source[x, y], {x, y} \[Element] mesh0, PlotPoints -> 75, that we will try and recover with the corresponding potential:

sol0 = NDSolveValue[{-Laplacian[u[x, y], {x, y}] ==
source[x, y] + NeumannValue[0, True] // Evaluate,
DirichletCondition[u[x, y] == 0, True]},  u, {x, y} \[Element] mesh0];
Plot3D[sol0[x, y], {x, y} \[Element] mesh0, PlotStyle -> Opacity[0.4],
PlotRange -> All, PlotTheme -> "ZMesh", PlotPoints -> 50] Let us sample this potential on a set of random points

data0 = RandomPoint[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]],500] // Sort;

ListPlot[data0, AspectRatio -> 1] and build the corresponding data set with the value of the potential on those points

data = Map[{#[], #[], sol0[#[], #[]]} &, data0];


Then the linear model follows from the knowledge of the data, data and the basis basis:

ff = Function[{x, y}, basis // Evaluate];
a = Map[ff @@ # &, Most /@ data];
a//Image//ImageAjust (looks a bit like the matrix) and we can fit the data as

fit[x_, y_] = basis.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)];


which is a pretty good fit!

Plot3D[fit[x, y] - sol0[x, y], {x, y} \[Element] mesh0,PlotRange -> All] Similarly we can invert for the source density

inv[x_, y_] =basis0.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)];
Plot3D[inv[x, y], {x, y} \[Element] mesh0, PlotRange -> All,
PlotTheme -> "ZMesh", PlotStyle -> Opacity[0.6]] Of course this inversion is a bit of an overkill to just get the density from the known potential BUT the framework works for any boundary condition and any sampling and arbitrary PDEs that mathematica can solve using FEM.

Indeed, compared to the analytic B-spline method, no extra work in needed to match the boundary conditions since the Mesh generator and FEM package takes care of that.

It is also worth pointing out that once a is known any data set can be inverted almost instantaneously.

To Do

1. I would be best to be able to define cubic splines as well on the mesh (following e.g. this).
2. One needs to write regularisation matrices on the mesh as well, in order to be able to invert ill-conditioned problems (but see this).
• @xzczd this is what I had in mind. Does it make sense to you? Mar 27, 2020 at 15:43
• You could use your code from this to speed up the NDSolve call here. You'd need to make a new coefficient that contains just the "LoadCoefficients" and re-use the stiffness matrix as that is the same in all cases. You possibly could even call LinrearSolve[stiffness, Method->"Pardiso"] once and apply that to each new right hand side. That will be very efficient. Keep going I am curious to see the result. It seems you are getting there. Mar 30, 2020 at 6:08
• @user21 thanks. I see your point in using the same stiffness matrix. Would you know how to get from the FEM framework the piece-wise linear basis (anchored on the vertices) which would be the analogue of basis0 above? I guess this is related to this question? mathematica.stackexchange.com/q/140601/1089 Mar 30, 2020 at 16:14
• @user21 I used to know this 15 years ago: I even did this in 4D and 6D (phase space) in order to solve for Boltzmann's equation on 4-simplices. But it was fairly inelegant and it would be much more useful to have an integrated solution like Mathematica's BSplineFunction Mar 30, 2020 at 17:24
• do you mean with shape functions like this? Mar 31, 2020 at 6:02

Thanks to @Henrik Schumacher's great help in extracting linear piecewise elements from FEM, let me provide a 1-spline solution appropriate for April's fool day. 2D case

reg = ImplicitRegion[(2 x^2 + y^2)^2 - 2 Sqrt x ( 2 x^2 - 3 y^2) + 2 (y^2 - x^2)<= 0, {x, y}]


and discretise it

R = ToElementMesh[R0=DiscretizeRegion[reg], MaxCellMeasure -> 0.015,
"MeshOrder" -> 1, MeshQualityGoal ->1]; R0 pts = R["Coordinates"]; n = Length[pts];
vd = NDSolveVariableData[
{"DependentVariables","Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {R}];
cdata = InitializePDECoefficients[vd, sd,"DiffusionCoefficients" ->
{{-IdentityMatrix}}, "MassCoefficients" -> {{1}}];
mdata = InitializePDEMethodData[vd, sd];


Discretisation yields

dpde = DiscretizePDE[cdata, mdata, sd];
stiffness = dpde["StiffnessMatrix"];
mass = dpde["MassMatrix"];


To see how it works, let us excite one basis element close to coordinate (0.4,0.1)

i = Nearest[pts -> "Index", {0.4, 0.1}][];
hatfun = ConstantArray[0., n];hatfun[[i]] = 1.;


This is how to interpolate it.

hatfuninterpolated = ElementMeshInterpolation[{R}, hatfun];
plot1 = Plot3D[hatfuninterpolated[x, y], {x, y} \[Element] R,
NormalsFunction -> None, PlotPoints -> 50, PlotTheme -> "Business",
BoxRatios -> {2, 1, 1}] In order to compute the corresponding potential let us extract the systemmatrix

bndplist =
Sort@DeleteDuplicates[Flatten[R["BoundaryElements"][[All, 1]]]];
intplist = Complement[Range[n], bndplist];


This is what DeployBoundaryConditions does to the stiffness matrix

systemmatrix = stiffness;
systemmatrix[[bndplist]] =
IdentityMatrix[n, SparseArray,
WorkingPrecision -> MachinePrecision][[bndplist]];


Factorizing the system matrix:

S = LinearSolve[systemmatrix, Method -> "Pardiso"];


Solving the actual equation yields the potential for this basis element.

solution = S[load];
solutioninterpolated = ElementMeshInterpolation[{R}, solution];
plot1 = Plot3D[solutioninterpolated[x, y] // Evaluate,
{x, y} \[Element] R,NormalsFunction -> None, PlotRange -> All,
ColorFunction ->
Function[{x, y, z}, RGBColor[1 - z/2, 1 - z, 1/2 - z]],
PlotTheme -> "Business", BoxRatios -> {2, 1, 1}] Let us now define a basis function

basis0 = Table[
hatfun = ConstantArray[0., n];
hatfun[[i]] = 1;
ElementMeshInterpolation[{R}, hatfun],
{i, 1, n}];


and compute its image

basis = Table[hatfun = ConstantArray[0., n];
ElementMeshInterpolation[{R}, solution],
{i, 1, n}];


Let us now pick a set of points on our fish

data0 = RandomPoint[R0, 1500] // Sort;
ListPlot[data0] and define a 'measured potential' from an (ad hoc random) set of 50 basis elements

hatfun0 = ConstantArray[0., n];
hatfun0[[RandomChoice[Range[n], 50]]] = 1;
sol0 = ElementMeshInterpolation[{R}, solution];
data = Map[{#[], #[], sol0[#[], #[]]} &, data0];


The linear model relating the basis to the data reads

ff = Function[{x, y}, Map[#[x, y] &, basis] // Evaluate];
a = Map[ff @@ # &, Most /@ data];

Clear[fit];
fit[x_, y_] := Module[{vec = Map[#[x, y] &, basis]},
vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];


Let us plot the fit:

Plot3D[fit[x, y] // Evaluate, {x, y} \[Element] R,
NormalsFunction -> None, PlotRange -> All,
ColorFunction ->
Function[{x, y, z}, RGBColor[1 - z/2, 1 - z, 1/2 - z]],
PlotTheme -> "Business", BoxRatios -> {2, 1, 1}] We can now also invert it:

Clear[inv];
inv[x_, y_] := Module[{vec = Map[#[x, y] &, basis0]},
vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];
Plot3D[inv[x, y] // Evaluate, {x, y} \[Element] R,
NormalsFunction -> None,
ColorFunction -> Function[{x, y, z},
RGBColor[1 - z/2, 1 - z, 1/2 - z]],
PlotTheme -> "Business", PlotPoints -> 50, BoxRatios -> {2, 1, 1},
PlotRange -> {0, 2}] It compares well with the input model:

hatfuninterpolated = ElementMeshInterpolation[{R}, hatfun0];
plot1 = Plot3D[hatfuninterpolated[x, y], {x, y} \[Element] R,
NormalsFunction -> None, PlotPoints -> 50, PlotTheme -> "Business",
BoxRatios -> {2, 1, 1},
PlotRange -> {0, 2}] Caveat: this is most likely not as efficient as it should be (see Henrik's comments). I could imagine e.g. that the way the basis function is defined is probably in part redundant w.r.t. to what is available within the FEM toolbox.

But it does illustrate that we can invert a given PDE with arbitrary sampling and ad hoc boundary condition on a set of linear piecewise basis function, which is differentiable, which is pretty cool IMHO. This question/answer provides means of regularising the inversion should this be needed (i.e. if a is poorly conditioned, with very small eigenvalues).

3D case

Let us give in one block the 3D code on a unit ball:

R = ToElementMesh[R0 = Ball[], MaxCellMeasure -> 0.125/16,
AccuracyGoal -> 1, "MeshOrder" -> 1];pts = R["Coordinates"];n = Length[pts];
vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u}, {x, y, z}}];
sd = NDSolveSolutionData[{"Space"} -> {R}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-IdentityMatrix}},
"MassCoefficients" -> {{1}}];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];
stiffness = dpde["StiffnessMatrix"];
mass = dpde["MassMatrix"];
bndplist = Sort@DeleteDuplicates[Flatten[R["BoundaryElements"][[All, 1]]]];
intplist = Complement[Range[n], bndplist]; systemmatrix = stiffness;
systemmatrix[[bndplist]] =
IdentityMatrix[n, SparseArray,
WorkingPrecision -> MachinePrecision][[bndplist]];
S = LinearSolve[systemmatrix, Method -> "Pardiso"];
basis0 = Table[
hatfun = ConstantArray[0., n];
hatfun[[i]] = 1;
ElementMeshInterpolation[{R}, hatfun],
{i, 1, n}];
basis = Table[
hatfun = ConstantArray[0., n];
hatfun[[i]] = 1; load = mass.hatfun;
solutioninterpolated = ElementMeshInterpolation[{R}, solution];
solutioninterpolated,
{i, 1, n}];

data0 = RandomPoint[R0, 500] // Sort;
hatfun0 = ConstantArray[0., n];
hatfun0[[RandomChoice[Range[n], 50]]] = 1;
sol0 = ElementMeshInterpolation[{R}, solution];

data = Map[{#[],#[],#[],sol0[#[], #[],#[]]} &, data0];
ff = Function[{x, y, z}, Map[#[x, y, z] &, basis] // Evaluate];
a = Map[ff @@ # &, Most /@ data];
Clear[fit];
fit[x_, y_, z_] := Module[{vec = Map[#[x, y, z] &, basis]},
vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];
Clear[inv];
inv[x_, y_, z_] := Module[{vec = Map[#[x, y, z] &, basis0]},
vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];


As a check let us look at the cross section through the mid-plane of the inverted density and the input density resp.

Plot3D[inv[x, y, 0] // Evaluate, {x, y} \[Element] Disk[],
NormalsFunction -> None, ColorFunction ->
Function[{x, y, z}, RGBColor[1 - z/2, 1 - z, 1/2 - z]],
PlotTheme -> "Business", PlotPoints -> 50, BoxRatios -> {1, 1, 1},
PlotRange -> {0, 2}]

hatfuninterpolated = ElementMeshInterpolation[{R}, hatfun0];
plot1 = Plot3D[hatfuninterpolated[x, y, 0], {x, y} \[Element] Disk[],
NormalsFunction -> None, PlotPoints -> 50, PlotTheme -> "Business",
BoxRatios -> {1, 1, 1},PlotRange -> {0, 2}] It seems to work fine!