# 1DPoisson equation with Dirac delta

I was writing my own FEM method to solve the Poisson equation

\begin{align} -u'' &= \exp(-c(x-1/2)^2)\\ u(0) &= u(1) = 0 \end{align}

where c=100.

and I'm reading from the book of Mats G. Larson, Fredrik Bengzon, "The Finite Element Method: Theory, Implementation, and Practice".

At page 42 there's the solution via an adaptive method of the equation. In my implementation, I obtain the same result, but the peak is higher, as in the following figure.

[

But, if I run my code with c=1000, I obtain the same solution of the picture of the book. So I suspect that the authors plotted the solution for c=1000 instead of c=100. I tried to solve the equation with Wolfram Mathematica, but I can't find it.

Can anyone confirm the solution for $$c=100$$?

This seems to confirm your suspicions that the plot in the book was obtained with $$c=1000$$:

ParametricNDSolveValue[{-u''[x] == Exp[-c (x - 1/2)^2], u[0] == 0, u[1] == 0}, u[x], {x, 0, 1}, c];
Plot[%[1000], {x, 0, 1}]


I also confirm that I obtained the plot you got for $$c=100$$.

• Many thanks for the answer, Marco! – VoB Nov 15 '19 at 21:58

MarcoB's answer is certainly the way to go to verify the result from the book. Nevertheless, I'd like to show a different approach that is geared towards verifying your FEM code. The finite element low level functions can help in doing that.

We start by setting up your PDE.

Needs["NDSolveFEM"]
sregion = Line[{{0}, {1}}];
vars = {x};
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u},
vars}];
sd = NDSolveSolutionData["Space" -> ToNumericalRegion[sregion]];

With[{c = 1000},
cdata = InitializePDECoefficients[vd, sd
, "DiffusionCoefficients" -> {{{{-1}}}}
, "LoadCoefficients" -> {{Exp[-c (x - 1/2)^2]}}];
]
bcdata = InitializeBoundaryConditions[vd,
sd, {{DirichletCondition[u[x] == 0, True]}}];


We use a first order mesh, as I assume that your code uses first order too.

mdata = InitializePDEMethodData[vd, sd,
Method -> {"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 1}}];


Here comes the benefit. DiscretizePDE can be used to both assemble the system matrices and to store the elements computed. We will look at those later.

(* save the computed elements and assemble the system matrices *)
ddata = DiscretizePDE[cdata, mdata, sd, "SaveFiniteElements" -> True,
"AssembleSystemMatrices" -> True];

dbcs = DiscretizeBoundaryConditions[bcdata, mdata, sd];

s = ddata["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbcs];
result = LinearSolve[s, l];


Now you can inspect the result from the LinearSolve step and generate an interpolating function.

if = ElementMeshInterpolation[{mdata["ElementMesh"]}, result];
Plot[if[x], {x, 0, 1}]


Looking at

ddata["LoadElements"]


and

ddata["StiffnessElements"]


Provides you with the computed values per element. You can use these to cross check your code. Or you could use NDSolveProcessEquations as an equation preprocessor for your own code:

femd = With[{c = 1000},
NDSolveProcessEquations[{-u''[x] == Exp[-c (x - 1/2)^2],
u[0] == 0, u[1] == 0}, u[x], {x, 0, 1},
Method -> "FiniteElement"][[1]]["FiniteElementData"]
];

femd["PDECoefficientData"]["DiffusionCoefficients"]
(* {{{{-1}}}} *)