# Normalization of solution of a PDE

I am trying to solve the following PDE:

pde = D[P[x, t], t] + D[J[x, t], x] == 0;

m1 = 0.6; m2 = 0.9;
J[x_, t_] := 1/2 ((-1 + 2 m1 - 2 x (-1 + m1 + m2)) P[x, t] + (-1 + x) x D[P[x, t], x]);
ic = P[x, 0] ==
Exp[-((x - 7./10)^2/(1./1000))]/(NIntegrate[Exp[-((x - 7/10)^2/(1/1000))], {x, 0.,1.}]); (*supposed to be a delta function but using an approximation*)
bc1 = (J[x, t] /. x -> 0) == 0.;
bc2 = (J[x, t] /. x -> 1) == 0.;
norm = Integrate[P[x, t], {x, 0, 1}] == 1;


Attempted solution:

1. DSolve failed to solve it. The analytical solution is here https://arxiv.org/pdf/q-bio/0508045.pdf (check page # 10 (iii))
2. NDSolve gave a solution. The approach used is as follows:
mol[n_Integer, o_ : "Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, "DifferenceOrder" -> o}}
molfem[measure_ : Automatic] := {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> MaxCellMeasure -> measure}};
nasol = NDSolveValue[{pde, ic, bc1, bc2}, P, {x, 0, 1}, {t, 0, 10}, Method -> mol[10000, 6]];
nasolfem = NDSolveValue[{pde, ic, bc1, bc2}, P, {x, 0, 1}, {t, 0, 10}, Method -> molfem[0.0001]];


But the solution isn't normalized: ListPlot[Table[{t, NIntegrate[nasol[x, t], {x, 0, 1}]}, {t, 0, 10, 1}]]

{{0, 1.}, {1, 0.993899}, {2, 0.971975}, {3, 0.946144}, {4, 0.920094}, {5, 0.894575}, {6, 0.869725}, {7, 0.845557}, {8, 0.82206}, {9, 0.799215}, {10, 0.777005}}


If we do the following (which is slow and needs lots of memory):

nasol = NDSolveValue[{pde, ic, bc1, bc2}, P, {x, 0, 1}, {t, 0, 10}, Method -> mol[1000000, 6]];
ListPlot[
Table[{t, NIntegrate[nasol[x, t], {x, 0, 1}]}, {t, 0, 10, 1}], PlotRange -> {0, 1}]


It's still not normalized.

Questions:

1. Is there an analytical solution?
2. How do I resolve the normalization issue?
• For #1: if you can assume that the dependence over x and t can be separated in the two functions, you could use the method of separation of variables: en.wikipedia.org/wiki/Separation_of_variables Jun 16, 2023 at 5:40

To answer this question we first reproduce solution from the paper Exact Solution of the Multi-Allelic Diffusion Model

pl = x^a (1 - x)^b JacobiP[l, a, b, 1 - 2 x];
cl = JacobiP[l, a, b, 1 - 2 x0] Gamma[2 R + l - 1] Gamma[
l + 1] (2 R + 2 l - 1)/Gamma[2 m1 + l]/Gamma[2 m2 + l];
a = 2 m1 - 1; b = 2 m2 - 1; R = m1 + m2; {m1 = 0.6, m2 = 0.9,
x0 = 0.7, n = 19};
P1 = Sum[cl pl Exp[-1/2 l (2 R + l - 1) t], {l, 0, n}];


Note that pl are eigenfunctions, while eigenvalues are given by lambda=1/2 l (2 R + l - 1).
Second, let's try to solve the Kolmogorov equation with Mathematica using DEigensystem as follows

n = 20;

op = With[{m1 = 6/10, m2 = 9/10},
1/2 D[((-1 + 2 m1 - 2 x (-1 + m1 + m2)) (p[x]) + (-1 + x) x D[p[x],
x]), x]];

{vals, funs} =
DEigensystem[{op, DirichletCondition[p[x] == 0, True]}, p[x], {x, 0, 1}, n];


We can compare vals with the analytical solution shown above as

With[{m1 = 6/10, m2 = 9/10}, Table[l (2 (m1 + m2) + l - 1)/2, {l, 0, n - 1}]]

(*Out[]= {0, 3/2, 4, 15/2, 12, 35/2, 24, 63/2, 40, 99/2, 60, 143/2, \
84, 195/2, 112, 255/2, 144, 323/2, 180, 399/2}*)

vals

(*Out[]= {0, 3/2, 4, 15/2, 12, 35/2, 24, 63/2, 40, 99/2, 60, 143/2, \
84, 195/2, 112, 255/2, 144, 323/2, 180, 399/2}*)


We see the perfect agreement. But eigenfunction funs are not normalized, so normalize it and compute coefficients c as follows

norm = Table[
NIntegrate[funs[[i]] funs[[i]], {x, 0, 1},

ff = Table[
NIntegrate[
funs[[i]]/Sqrt[norm[[i]]] funs[[j]]/Sqrt[norm[[j]]], {x, 0, 1},
Method -> "LocalAdaptive"], {i, n}, {j, n}];

dd = Integrate[DiracDelta[x - 7/10] funs, {x, 0, 1}] // N; dn =
Table[dd[[i]]/Sqrt[norm[[i]]], {i, Length[dd]}];

c = LinearSolve[ff, dn];


Finally, we define the analytical solution in a form

P = Sum[c[[i]] funs[[i]]/Sqrt[norm[[i]]] Exp[-vals[[i]] t], {i,
n}];
normP = NIntegrate[P /. t -> 0, {x, 0, 1}, Method -> "LocalAdaptive"];


Now we can compare two analytical solutions P (blue solid line), and P1(red dashed line) in one plot for different times as

Table[Plot[Evaluate[{P/normP, P1} /. t -> t0], {x, 0, 1},
PlotLabel -> t0, PlotRange -> All,
PlotStyle -> {{Blue}, {Red, Dashed}}], {t0, .0, 1.5, .1}]


We also can plot these data in 3D

{Plot3D[P, {x, 0, 1}, {t, 0.05, 1.5}, Mesh -> None,
ColorFunction -> Hue, PlotRange -> All],
Plot3D[P1, {x, 0, 1}, {t, 0.05, 1.5}, Mesh -> None,
ColorFunction -> Hue, PlotRange -> All]}


Check normalization

Table[NIntegrate[P1, {x, 0, 1}, Method -> "LocalAdaptive"], {t, 0, 10,
1}] // Quiet
{1.00007, 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}

Table[NIntegrate[P/normP, {x, 0, 1}, Method -> "LocalAdaptive"], {t,
0, 10, 1}] // Quiet

(*Out[]= {1., 0.999981, 0.999981, 0.999981, 0.999981, 0.999981, \
0.999981, 0.999981, 0.999981, 0.999981, 0.999981}*)


Update 1 Let's consider the numerical solution to this problem with Mathematica FEM. First, we prepare handmade mesh as follows

Needs["NDSolveFEM"]

xmin = 0.; xmax = 1.; reg =
ImplicitRegion[xmin <= x <= xmax, {x}];
mesh =ToElementMesh[reg, MaxCellMeasure -> 1. 10^-3];
mesh1 = ToElementMesh[reg, MaxCellMeasure -> 1. 10^-4];
mesh2 = ToElementMesh[reg, MaxCellMeasure -> 1. 10^-5];


Next we compute 3 solutions with using mesh,mesh1,mesh2 and initial condition p[x, 0] == P /. t -> .05, we have

pde = D[p[x, t], t] + D[J[x, t], x] == 0;
J[x_, t_] :=
1/2 ((-1 + 2 m1 - 2 x (-1 + m1 + m2)) p[x, t] + (-1 + x) x D[p[x, t],
x]);
ic = p[x, 0] == P /. t -> .05; sol =
NDSolve[{pde, ic,
DirichletCondition[p[x, t] == 0, x == xmin || x == xmax]}, p,
Element[{x}, mesh], {t, 0, 1.5}];sol1 =
NDSolve[{pde, ic,
DirichletCondition[p[x, t] == 0, x == xmin || x == xmax]}, p,
Element[{x}, mesh1], {t, 0, 1.5}];sol2 =
NDSolve[{pde, ic,
DirichletCondition[p[x, t] == 0, x == xmin || x == xmax]}, p,
Element[{x}, mesh2], {t, 0, 1.5}];


Finally we show in one plot analytical solution P(1), and numerical solutions sol,sol1,sol2- 2,3,4 as

Plot[Evaluate[{P,
p[x, t - .05] /. {sol[[1]], sol1[[1]], sol2[[1]]}} /.
t -> 1.5], {x, 0, 1}, PlotLegends -> Automatic]


We can see from this picture, that the numerical solution approaches the analytical solution with an increasing number of mesh elements. Also, we can compute the norm for every solution at t=1.45

normp =
NIntegrate[p[x, 1.45] /. sol[[1]], {x, 0, 1},

(*Out[]= 0.988646*)
normp1 =
NIntegrate[p[x, 1.45] /. sol1[[1]], {x, 0, 1},

(*Out[]= 0.99864*)
normp2 =
NIntegrate[p[x, 1.45] /. sol2[[1]], {x, 0, 1},

(*Out[]= 1.00392*)


Update 2. We also can compute eigenfunctions and eigenvalues with using NDEigensystem as follows

n=20;
op1 = With[{m1 = 6/10,
m2 = 9/10}, {-1/
2 D[((-1 + 2 m1 - 2 x (-1 + m1 + m2)) (p[x]) + (-1 + x) x D[p[x],
x]), x], DirichletCondition[p[x] == 0, True]}]; {vals1,
funs1} =
NDEigensystem[op1, p[x], {x, 0, 1}, n,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}},
"Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];


Using these functions we can define numerical solutions in a form (code is same as above for analytical solution)

norm = Table[
NIntegrate[funs1[[i]]^2, {x, 0, 1}, Method -> "LocalAdaptive"], {i,
Length[vals1]}]

ff =
Table[NIntegrate[
funs1[[i]]/Sqrt[norm[[i]]] funs1[[j]]/Sqrt[norm[[j]]], {x, 0, 1},
Method -> "LocalAdaptive"], {i, n}, {j, n}];
dd = Integrate[DiracDelta[x - 7/10] funs1, {x, 0, 1}] // N; dn =
Table[dd[[i]]/Sqrt[norm[[i]]], {i, Length[dd]}];
c = LinearSolve[ff, dn];
P = Sum[c[[i]] funs1[[i]]/Sqrt[norm[[i]]] Exp[vals[[i]] t], {i, n}];


Visualization

Plot3D[Evaluate[P], {x, 0, 1}, {t, 0.05, 1.5}, Mesh -> None,
ColorFunction -> Hue, PlotRange -> All, PlotPoints -> 100]


This picture looks very similar to the analytical solution, but if we compute the norm, then we have large discrepancies with the analytical solution,

Table[NIntegrate[P, {x, 0, 1}, Method -> "LocalAdaptive"] //
Quiet, {t, 0, 10, 1}]

(*{1.01157, 1.00819, 0.997508, 0.984611, 0.971379, 0.958216, 0.945208, 0.932372, 0.919709, 0.907218, 0.894896}*)


The problem is that vals1 computed with a high error compare to the analytical one,

vals1

(* {-0.0136747, -1.53703, -4.06629, -7.60031, -12.1385, -17.6803, -24.2255, -31.774, -40.3254, -49.8796, -60.4365, -71.9961, -84.5582, -98.1227, -112.69, -128.259, -144.83, -162.404, -180.98, -200.557}*)
-vals
{0, -(3/2), -4, -(15/2), -12, -(35/2), -24, -(63/2), -40, -(99/
2), -60, -(143/2), -84, -(195/2), -112, -(255/2), -144, -(323/
2), -180, -(399/2)}

• Thanks a lot! I really appreciate it. I actually have lots of questions but let me start with two basic ones: 1. How were you able to read the eigenfunction and eigenvalue from the solution? 2. Could you please explain how you were able to formulate the PDE problem as Eigensystem problem?
– sra
Jun 19, 2023 at 5:19
• @sra Please, pay attention that in the paper you cited they expressed solution in a form $P=\sum_{l=0}^{\infty}c_l p_l(x)e^{-\lambda_l t}$. If we put this expression in PDE, then we come to the eigenvalues problem as it mentioned above. Your first question is not so clear. Do you interested in how to separate $p_l, \lambda_l$ from analytical solution? Jun 19, 2023 at 5:39
• Thanks a lot. I think I understand eigenvalue formulation now. However, I still do not get DEigensystem. Let's say that if I change m1 = m2=0.5 then mathematica fails to solve it. It has to do with the DirichletCondition because for m1-m2 p[x] != 0 at x=0 and x=1. Can DEigensystem be used to solve it for general m1 and m2?
– sra
Jun 19, 2023 at 18:02
• @sra No, DEigensystem can't be used for arbitrary m1,m2. It is why I add section Update 2. Jun 20, 2023 at 1:42
• For m1 = m2  the p[x]==0 boundary condition (BC) does not work. For example: op = With[{m1 = 5/10, m2 = 5/10}, 1/2 D[((-1 + 2 m1 - 2 x (-1 + m1 + m2)) (p[x]) + (-1 + x) x D[p[x],x]), x]]; {vals, funs} = DEigensystem[{op, DirichletCondition[p[x] == 0 , True]}, p[x], {x, 0, 1}, n]; How do we impose BC in such case?
– sra
Jun 20, 2023 at 3:44

molfem[measure_ : Automatic] := {"MethodOfLines","SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> MaxCellMeasure -> measure}};

nasolfem =NDSolveValue[{pde, ic, bc1, bc2}, P, {x, 0, 1}, {t, 0, 10} ,Method-> molfem[0.01 ] ];
Plot3D[nasolfem[x, t], {x, 0, 1}, {t, 0, 10}, PlotRange -> All,AxesLabel -> {"x", "t", "P[x,t]"}]


scaling function

mw =.
mw[t_?NumericQ] :=Block[{x},NIntegrate[nasolfem[x, t], {x, 0, 1}, Method -> "FiniteElement"]]
Plot[mw[t], {t, 0, 10}, PlotLabel -> "Scaling function mw[t]"]


normalized solution

Show[Table[
ParametricPlot3D[Evaluate[{x, t, nasolfem[x, t]/mw[t]}], {x, 0, 1} ,BoxRatios -> {1, 1, 1}, PlotStyle -> Hue[t ]] , {t,Subdivide[0, 1.5, 15 3] // Rest}], PlotRange -> {{0, 1}, {0, 1.5}, {0, 5}},
PlotLabel ->"normalized solution \!$$\*FractionBox[\(P[x, t]$$,$$mw[t]$$]\)", AxesLabel -> {x, t, P/mw}]


• Thank you very much (+1). :) Jun 17, 2023 at 11:12
• @AlexTrounev You 're welcome. FEM solution only roughly approximates your exact solution. Why these differences? Jun 17, 2023 at 11:32
• Use handmade mesh and DirichletCondition[P[x, t] == 0, True]. Jun 17, 2023 at 12:47
• @AlexTrounev Thanks, that makes it easier Jun 17, 2023 at 14:06
• @UlrichNeumann Thanks a lot for helping out.
– sra
Jun 19, 2023 at 5:18