1
$\begingroup$

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}]

enter image description here

It's still not normalized.

Questions:

  1. Is there an analytical solution?
  2. How do I resolve the normalization issue?
$\endgroup$
1

2 Answers 2

2
$\begingroup$

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}, 
    Method -> "LocalAdaptive"], {i, n}];

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}]

Figure 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]}

Figure 2

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]

Figure 3

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}, 
  Method -> "LocalAdaptive"]

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

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

(*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]

Figure 4

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)}
$\endgroup$
10
  • $\begingroup$ 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? $\endgroup$
    – sra
    Jun 19, 2023 at 5:19
  • 1
    $\begingroup$ @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? $\endgroup$ Jun 19, 2023 at 5:39
  • $\begingroup$ 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? $\endgroup$
    – sra
    Jun 19, 2023 at 18:02
  • $\begingroup$ @sra No, DEigensystem can't be used for arbitrary m1,m2. It is why I add section Update 2. $\endgroup$ Jun 20, 2023 at 1:42
  • $\begingroup$ 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? $\endgroup$
    – sra
    Jun 20, 2023 at 3:44
2
$\begingroup$

answer #2(only FEM-solution works!)

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]"}]

enter image description here

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]"]

enter image description here

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}]

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.