# Solving the Lotka-McKendrick model with NDSolve

The Lotka-McKendrick model is a demographic model that represents the way a population changes over time due to fertility and mortality. For an age-specific population density $$u(a, t)$$, and a total birth rate $$\Lambda(t)$$, the following equations should be satisfied:

$$\begin{eqnarray} \frac{\partial u}{\partial a} + \frac{\partial u}{\partial t} & = & -\mu(a) u(a, t) \\ \Lambda(t) & = & u(0, t) = \int_{0}^{\infty} da\,u(a,t) f(a,t) \\ u(a, 0) & = & u_0(a) \end{eqnarray}$$

Here, $$\mu(a)$$ is an age-specific force of mortality, $$f(a)$$ is an age-specific fertility rate, and $$u_0(a)$$ is an initial condition.

If it weren't for the integral in the boundary condition $$\Lambda(T) = u(0, t)$$, we'd be home free. In fact, DSolve would suffice, using the standard technique of integrating along the characteristic lines of the first-order PDE:

lkPDE = {D[u[a, t], a] + D[u[a, t], t] == -\[Mu][a]*u[a, t],
u[a, 0] == u0[a]};

DSolve[lkPDE, u[a, t], {a, t}]
(* {{u[a, t] ->
E^(Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a}] - Inactive[Integrate][-\[Mu][K[1]],
{K[1], 1, a - t}])*u0[a - t]}} *)



Sticking the integral in there makes everything fall apart.

lkIntegral =
u[0, t] == Integrate[u[x, t]*f[x], {x, 0, Infinity}];

DSolve[Flatten@{lkPDE, lkIntegral}, u[a, t], {a, t}]
(* returns unevaluated *)


You can write an analytic solution down, but as an alternative, I would like to use NDSolve, especially as numerical methods will generalize to cases where the analytic solutions don't exist or are too complicated to be useful.

Sadly, NDSolve also chokes, even with appropriate concessions to reality. Let's choose very simple parameters:

$$\begin{eqnarray} \mu(a) & = & 1/80 \\ f(a) & = & \left\{ \begin{array} & 1/10 & 20 \le a < 45 \\ 0 & \text{otherwise} \end{array} \right. \end{eqnarray}$$

Even so, we need a simpler integral condition because Integrate isn't quite smart to handle that piecewise function.

simpleLkIntegral =
u[0, t] == Integrate[u[x, t], {x, 20, 45}]/10

NDSolve[{
lkPDE /. \[Mu][_] -> 1/80 /. u0[a_] :> 1/80,
simpleLkIntegral
},
u,
{a, 0, 100},
{t, 0, 100}]
(* returns unevaluated, with the an NDSolve::litarg message complaining about the integral *)


In order to appease NDSolve::litarg, I try rewriting the integral with a replacing x as the variable of integration, which yields no joy:

simpleLkIntegral2 =
u[0, t] == Integrate[u[a, t], {a, 20, 45}]/10

NDSolve[{
lkPDE /. \[Mu][_] -> 1/80 /. u0[a_] :> 1/80,
simpleLkIntegral2
},
u,
{a, 0, 100},
{t, 0, 100}]
(* returns unevaluated, with a
NDSolve::overdet: There are fewer dependent variables, {u[a,t]}, than equations, so the system is overdetermined.
*)


At this point, I feel like I've more or less run out of road, but was wondering if there was some way to force NDSolve to do what I want. UPDATE: I did try the model again with a different set of initial conditions, ones which allow for consistency between the boundary and initial conditions from $$t = 0$$ on, as shown below:

simpleLkInit = With[{m = 1/80},
u0[a_] :> Piecewise[{{c - m*a, 0 <= a <= c/m}}, 0]];

simpleLkNormalization = First@Solve[
{simpleLkIntegral2 /. t -> 0 /. u[a_, 0] :> u0[a] /. simpleLkInit,
c > 0}, c]
(* c -> 65/96 *)


Plugging this into NDSolve gives the same issue with overdetermination (so presumably the consistency of the boundary condition is never even checked):

NDSolve[{lkPDE /. \[Mu][_] -> 1/80 /. simpleLkInit /.
simpleLkNormalization, simpleLkIntegral2}, u, {a, 0, 100}, {t, 0,
100}]
(* Unevaluated, with NDSolve::overdet message *)


While the strategy of discretizing the system in age manually, as in Chris K's fine answer, is entirely viable, this essentially boils down to using the method of lines, which is the approach NDSolve itself uses. I would like to see if NDSolve itself can do the discretization, or at least if I can use it to drive move of the problem.

• We can solve it, but with message: NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. Commented Sep 19, 2020 at 10:52
• So, the solution is supposed to be unsmooth somewhere? (The initial condition is enough to determine the classic solution. ) Commented Sep 20, 2020 at 14:01
• @xzczd It is not only an unsmooth solution but not unique as well. Commented Sep 20, 2020 at 15:45
• @AlexTrounev The question still needs clarification in my view. Since we're not looking for a classic solution, additional condition is needed to determine a particular solution, then what's it? I guess it's certain continuity or conservation. (Just like the heat continuity condition for heat conduction equation. ) But I can't go further given that I'm not familiar with the underlying model. Commented Sep 21, 2020 at 2:20
• @xzczd It is also known as the von Foerster equation. It's a standard model in mathematical demography. It includes three processes: aging ($\partial u/\partial a$), death ($-\mu(a)u$), and reproduction ($\int f(a)u(a) da$). Despite the weird boundary condition, it is linear. You can read more about it here. Commented Sep 23, 2020 at 17:35

I'm not an expert on age-structured populations (particularly this continuous-time model) and I know better numerical methods exist, but why not just discretize in age a and solve the resulting large system of ODEs?

(NB: doublecheck the details of my discretization if you use this for anything serious; I wasn't too careful in how I put the da's in!)

imax = 100;
da = 1.0;

f[a_] := If[20 <= a < 45, 1/10, 0];
μ[a_] := 1/80;
u0[a_] := If[a <= 10, 1/80, 0];

eqns = Join[
{u[0]'[t] == -μ[0] u[0][t] - u[0][t]/da + Sum[u[i][t] f[i da], {i, 0, imax}]},
Table[u[i]'[t] == -(u[i][t] - u[i - 1][t])/da - μ[i da] u[i][t], {i, 1, imax}]
];
ics = Table[u[i][0] == u0[i da], {i, 0, imax}];
unks = Table[u[i], {i, 0, imax}];

tmax = 160;
sol = NDSolve[{eqns, ics}, unks, {t, 0, tmax}][[1]];

frames = Table[
ListPlot[Table[{i da, u[i][t]}, {i, 0, imax}] /. sol,
PlotRange -> {0, 0.06}, PlotLabel -> t, AxesLabel -> {"a", "u"}]
, {t, 0, tmax}];

ListAnimate[frames]


I started with a step-function of u0[a] to illustrate a few things:

1. You can see the population distribution move to the right as individuals age.
2. There's a baby boom when the initial population goes through reproductive ages 20-45, and also echoes as their kids reproduce.
3. The population approaches a stable age distribution, then grows exponentially.
4. Somewhat problematic: the discretization of the advection term results in numerical diffusion, blurring the initial step-function distribution over time. Higher resolution (smaller da) helps, and if you're only interested in the long-term behavior or smooth age-distributions, I think this isn't too bad. This is where more sophisticated numerical methods may help.

Finally, an advantage of this approach is you can look at the eigenvalues and eigenvectors to get more info. Linearizing to make a matrix A:

A = D[eqns /. ((_ == rhs_) -> rhs) /. (var_[t] -> var), {unks}];
{λ, v} = Eigensystem[A];
λ[[-3 ;; -1]]

(* {-0.0370978 + 0.184096 I, -0.0370978 - 0.184096 I, 0.0163063 + 0. I} *)


The last eigenvalue is the dominant, which gives the asymptotic growth rate to be 0.0163063 per year. The subdominant eigenvalues are complex; I think the imaginary part gives the approximate period of those baby boom echoes:

Abs[2 π/Im[λ[[-2]]]]
(* 34.1299 *)


Finally, the eigenvector associated with the dominant eigenvalue gives the stable age distribution:

ListPlot[Abs[v[[-1]]]]


EDIT:

Better yet, you can just convert this to a discrete-time, discrete-state Leslie matrix model. As long as the time step matches the size of the age classes, there is no spurious numerical diffusion.

Make the Leslie matrix:

L = SparseArray[
Join[
Table[{1, i + 1} -> f[i da] da, {i, 0, imax}],
Table[{i + 2, i + 1} -> 1 - μ[i da] da, {i, 0, imax - 1}]
], {imax + 1, imax + 1}
];


Project forward in time:

n = Table[If[i <= 11, 1/80, 0], {i, imax + 1}];
res = Join[{n}, Table[n = L.n, {t, 1, tmax}]];

frames = Table[
ListPlot[res[[t + 1]], PlotLabel -> t da, PlotRange -> {0, da 0.06}, AxesLabel -> {"a", "u"}]
, {t, 0, tmax/da}];
ListAnimate[frames]


Asymptotic growth rate checks out:

Log[Max[Re[Eigenvalues[A]]]]/da
(* 0.0162194 *)


EDIT 2:

I think you'll end up stuck with manual discretization in age, because the boundary condition is so weird compared to most typical PDEs. The discrete time-step in my matrix approach avoids numerical diffusion, which is important to maintain the shape if there are steps in the initial conditions (this should be a stringent test for any answer that tries to solve this problem).

The only thing I have left to offer is to force NDSolve to solve the continuous-time system in the same way as the discrete-time version using Method->"ExplicitEuler" and step size equal to the width of an age class. (note I had to tweak my discretization a little bit).

Here's a nice hi-resolution example:

imax = 1000;
da = 0.1;

f[a_] := If[20 <= a < 45, 1/10, 0];
μ[a_] := 1/80;
u0[a_] := If[a < 5, 0.1, 0];

eqns = Join[
{u[0]'[t] == -μ[0] u[0][t] - u[0][t]/da + Sum[u[i][t] f[i da], {i, 0, imax}]},
Table[u[i]'[t] == -(u[i][t] - u[i - 1][t])/da - μ[(i - 1) da] u[i - 1][t], {i, 1, imax}]
];
ics = Table[u[i][0] == u0[i da], {i, 0, imax}];
unks = Table[u[i], {i, 0, imax}];

tmax = 160;
sol = NDSolve[{eqns, ics}, unks, {t, 0, tmax},
Method -> "ExplicitEuler", StartingStepSize -> da][[1]];

frames = Table[
ListPlot[Table[{i da, u[i][t]}, {i, 0, imax}] /. sol,
PlotRange -> {0, 0.2}, PlotLabel -> t, AxesLabel -> {"a", "u"},
Joined -> True]
, {t, 0, tmax}];
ListAnimate[frames]


• Nice attempt (+1)! Your solution corresponds to my first solution with the message NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. Commented Sep 19, 2020 at 20:01
• @AlexTrounev Is that warning really something to worry about? I mean, the dynamics of the system seems to take care of it naturally. Commented Sep 19, 2020 at 20:19
• Yes, it looks natural, but NDSolve does not agree with this. :) Commented Sep 19, 2020 at 21:30
• @AlexTrounev In any case, I don't get that error with my code! Commented Sep 20, 2020 at 2:04
• You don't get error since you solved system of ODEs, not PDE. Commented Sep 20, 2020 at 10:57

There is no unique solution for data provided by @Pillsy, since boundary and initial conditions are inconsistent. To show it we just use exact solution in a form:

\[Mu][a_] := 1/80; u0[a_] := 1/80;
u[a_, t_] :=
E^(Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a}] -
Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a - t}])*u0[a - t]
u[0, t_] := Integrate[u[x, t], {x, 20, 45}]/25;


Now we can plot u[0,t] as follows:

Plot[u[0, t], {t, 0, 100}]


So it is smooth function and we can make interpolation in a form

lst = Table[{t, u[0, t] // N}, {t, 0, 100, 1}];

ut = Interpolation[lst];


With ut we can use NDSolve directly

sol = NDSolveValue[{D[v[a, t], a] + D[v[a, t], t] == -\[Mu][a]*
v[a, t], v[a, 0] == u0[a], v[0, t] == ut[t]},
v, {a, 0, 100}, {t, 0, 100}]


Here we got message NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. Nevertheless numerical solution can be plot and it looks like waved function

Plot3D[sol[a, t], {a, 0, 100}, {t, 0, 100}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic]


To avoid the inconsistent of boundary and initial conditions we put in the beginning of the code

u[0, t_] := Integrate[u[x, t], {x, 20, 45}]/25;


Then we get smooth numerical solution

Now we can use method of line implemented by Chris K with some appropriate modifications

Clear[u];

imax = 200;
da = 1/2;

f[a_] := If[20 <= a < 45, 1/25, 0];
\[Mu][a_] := 1/80;
u0[a_] := 1/80;

eqns = Join[{u[0]'[t] ==
da/2 Sum[(u[i + 1]'[t] f[(i + 1) da] + u[i]'[t] f[i da]), {i, 0,
imax - 1}]},
Table[u[i]'[
t] == -(u[i][t] - u[i - 1][t])/da - \[Mu][i da] u[i][t], {i, 1,
imax}]];
ics = Table[u[i][0] == u0[i da], {i, 0, imax}];
unks = Table[u[i], {i, 0, imax}];

tmax = 100;
sol1 = NDSolve[{eqns, ics}, unks, {t, 0, tmax}][[1]];


We make visualization of numerical solution of the system of ODEs and see identical picture as we got for PDE

ListPlot3D[
Flatten[Table[{i da, t, Evaluate[u[i][t] /. sol1]}, {i, 0, imax}, {t,
0, tmax, 1}], 1], AxesLabel -> {"a", "t", "u"},
ColorFunction -> "Rainbow", PlotRange -> All, Mesh -> None]


And more detailed comparison also shows agreement of two solutions

Show[Plot[{sol[10, t]}, {t, 0, 100}],
ListPlot[Table[{t, u[20][t] /. sol1}, {t, 0, 100, 2}],
PlotStyle -> Red]]


One example of growing population with consistent boundary and initial conditions:

\[Mu][a_] := 1/80; u0[a_] := Exp[-a/45];
f[a_] := Piecewise[{{1/10/1.2298542626633067, 20 <= x < 45}, {0,
True}}];
ue[a_, t_] :=
E^(Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a}] -
Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a - t}])*u0[a - t]
u1[t_] := NIntegrate[ue[x, t] f[x], {x, 0, 100}] // Quiet;

lst = Table[{t, u1[t]}, {t, 0, 100, 1}];

ut = Interpolation[lst];

sol = NDSolveValue[{D[v[a, t], a] + D[v[a, t], t] == -\[Mu][a]*
v[a, t], v[a, 0] == u0[a], v[0, t] == ut[t]},
v, {a, 0, 100}, {t, 0, 100}]


Visualisation

Plot3D[sol[a, t], {a, 0, 100}, {t, 0, 100}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic]


And same solution with method of lines:

imax = 500;
da = 100/imax;

f[a_] := If[20 <= a < 45, 1/10/1.2298542626633067, 0];
\[Mu][a_] := 1/80;
u0[a_] := Exp[-a/45];

eqns = Join[{u[0]'[t] ==
da/2 Sum[(u[i + 1]'[t] f[(i + 1) da] + u[i]'[t] f[i da]), {i, 0,
imax - 1}]},
Table[u[i]'[
t] == -(u[i][t] - u[i - 1][t])/da - \[Mu][i da] u[i][t], {i, 1,
imax}]];
ics = Table[u[i][0] == u0[i da], {i, 0, imax}];
unks = Table[u[i], {i, 0, imax}];

tmax = 100;
sol1 = NDSolve[{eqns, ics}, unks, {t, 0, tmax}][[1]];


Let compare two solution and find out that they have small discrepancies (due to large da)

Table[Show[
Plot[{sol[i da, t]}, {t, 0, 100}, AxesLabel -> Automatic,
PlotLabel -> Row[{"a = ", i da}]],
ListPlot[Table[{t, u[i][t] /. sol1}, {t, 0, 100, 2}],
PlotStyle -> Red]], {i, 0, imax, 20}]


Last example provided by Pillsy shows divergence of two solutions found out with two methods even if initial data and boundary conditions are consistent. First method:

\[Mu][a_] := 1/80; u0[a_] := If[0 <= a <= 325/6, 65/96 - a/80, 0];
f[a_] := Piecewise[{{1/10, 20 <= x < 45}, {0, True}}];
ue[a_, t_] :=
E^(Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a}] -
Inactive[Integrate][-\[Mu][K[1]], {K[1], 1, a - t}])*u0[a - t]
u1[t_] := NIntegrate[ue[x, t], {x, 20, 45}]/10 // Quiet;

lst = Table[{t, u1[t]}, {t, 0, 100, 1/4}];

ut = Interpolation[lst];

sol = NDSolveValue[{D[v[a, t], a] + D[v[a, t], t] == -\[Mu][a]*
v[a, t], v[a, 0] == u0[a], v[0, t] == ut[t]},
v, {a, 0, 100}, {t, 0, 100}];


Second method:

imax = 500;
da = 100/imax;

f[a_] := If[20 <= a < 45, 1/10, 0];
\[Mu][a_] := 1/80;
u0[a_] := If[0 <= a <= 325/6, 65/96 - a/80, 0];

eqns = Join[{u[0]'[t] ==
da/2 Sum[(u[i + 1]'[t] f[(i + 1) da] + u[i]'[t] f[i da]), {i, 0,
imax - 1}]},
Table[u[i]'[
t] == -(u[i][t] - u[i - 1][t])/da - \[Mu][i da] u[i][t], {i, 1,
imax}]];
ics = Table[u[i][0] == u0[i da], {i, 0, imax}];
unks = Table[u[i], {i, 0, imax}];

tmax = 100;
sol1 = NDSolve[{eqns, ics}, unks, {t, 0, tmax},
Method -> {"EquationSimplification" -> "Residual"}][[1]];


Now we plot solutions together and see divergence

Table[Show[

Plot[{sol[i da, t]}, {t, 0, 100}, AxesLabel -> Automatic,
PlotLabel -> Row[{"a = ", i da}]],
ListPlot[Table[{t, u[i][t] /. sol1}, {t, 0, 100, 2}],
PlotStyle -> Red, PlotRange -> All]], {i, 0, imax, 20}]


Nevertheless, we can consider all tests above as verification of numerical method of lines for this problem. Now we make next step to develop code with known error of $$h^4$$, where $$h$$ is step size. For this we use function GaussianQuadratureWeights[] to integrate with n-point Gaussian formula for quadrature, and function FiniteDifferenceDerivative for approximation of derivative $$\frac {\partial u}{\partial x}$$ with DifferenceOrder of 4. First we call utilities:

Needs["DifferentialEquationsNDSolveProblems"]
Needs["DifferentialEquationsNDSolveUtilities"]
Get["NumericalDifferentialEquationAnalysis"]


Second step, we define derivative matrix m and integral vector int:

np = 400; g = GaussianQuadratureWeights[np, 0, 100];
ugrid = g[[All, 1]]; weights = g[[All, 2]];
fd = NDSolveFiniteDifferenceDerivative[Derivative[1], ugrid]; m =
fd["DifferentiationMatrix"]; vart =
Table[u[i][t], {i, Length[ugrid]}]; vart1 =
Table[u[i]'[t], {i, Length[ugrid]}]; ux = m.vart; var =
Table[u[i], {i, Length[ugrid]}];
f[a_] := If[20 <= a < 45, 1/10, 0]; int =
Table[f[ugrid[[i]]] weights[[i]], {i, np}];
\[Mu][a_] := 1/80;
u0[a_] := If[0 <= a <= 325/6, 65/96 - a/80, 0];


Third step, we define the system of equations:

eqns = Join[{D[u[1][t], t] == int.vart1},
Table[D[u[i][t], t] == -ux[[i]] - \[Mu][ugrid[[i]]] u[i][t], {i, 2,
Length[ugrid]}]];
ics = Table[u[i][0] == u0[ugrid[[i]]], {i, Length[ugrid]}];


Finally we solve system as

tmax = 100;
sol1 = NDSolve[{eqns, ics}, var, {t, 0, tmax},
Method -> {"EquationSimplification" -> "Residual"}];


With this code we made research to check how solution diverges with np increasing:

{np, {u[1][100] /. sol1[[1]], u[np][100] /. sol1[[1]]}}

{100, {4.0455, 0.197089}}

{200, {3.791317314610565, 0.19572819660924937}};

{400, {3.6951293716506926, 0.1949809561721866}};

{800, {3.70082201902361, 0.19456320959442788}};



So numerical solution converges very fast with np increasing, and for np=800 we have picture

lst1 = Flatten[
Table[{t, ugrid[[i]], u[i][t] /. sol1[[1]]}, {t, 0, 100, 2}, {i, 1,
Length[ugrid], 5}], 1];
ListPlot3D[lst1, Mesh -> None, ColorFunction -> Hue, PlotRange -> All,
AxesLabel -> {"t", "a"}]


We have run several tests with known exact solution and got a good agreement of the exact and numerical solution getting with the last code. Example 1 from Numerical methods for the Lotka–McKendrick’s equation (there are typos in this paper in equations (6.8),(6,9) I have corrected using Mathematica 12.1):

f[a_]:=2; \[Mu][a_] := 1/(1 - a);
p0[x_] := If[x <= 1/2, (1 - 2 x)^3 (1 - x), 31 (2 x - 1)^3 (1 - x)];
u0[a_] := p0[a];


With this data we have to compute u[a,t] on {a,0,1},{t,0,1} and then compare exact solution $$u(0,t)=B(t)$$ with numerical solution:

B[t_] := If[t <= 1/2,
217 - 186 t - 372 t^2 - 248 t^3 - 216 E^t Cos[t] + 396 E^t Sin[t],
1/(Sqrt[E] (Cos[1/2]^2 + Sin[1/2]^2)) (-7 Sqrt[E] Cos[1/2]^2 +
6 Sqrt[E] t Cos[1/2]^2 + 12 Sqrt[E] t^2 Cos[1/2]^2 +
8 Sqrt[E] t^3 Cos[1/2]^2 - 216 E^(1/2 + t) Cos[1/2]^2 Cos[t] +
768 E^t Cos[t] Sin[1/2] - 7 Sqrt[E] Sin[1/2]^2 +
6 Sqrt[E] t Sin[1/2]^2 + 12 Sqrt[E] t^2 Sin[1/2]^2 +
8 Sqrt[E] t^3 Sin[1/2]^2 - 216 E^(1/2 + t) Cos[t] Sin[1/2]^2 -
768 E^t Cos[1/2] Sin[t] + 396 E^(1/2 + t) Cos[1/2]^2 Sin[t] +
396 E^(1/2 + t) Sin[1/2]^2 Sin[t])];


In Figure 10 shown numerical solution (left) and exact solution (right, blue line) with numerical solution (red points): Finally we have to compare absolute error $$|B(t)-u(0,t)|$$ for code by Chris K and my code to find out accuracy of two code. For Chris code it is obvious that error is of h and for my code theoretically it should be of $$h^3$$. But since we solve PDE it is not so perfect. In Figure 11 shown absolute error for Chris code (left) and for my code (right) for imax=np=800. It looks like my code has error of $$h^{3/2}$$ not $$h^3$$. To improve the accuracy of my code we add point a=0 to the grid for differentiation matrix and finally have

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"]; \
Get["NumericalDifferentialEquationAnalysis"];

np = 40; g = GaussianQuadratureWeights[np, 0, 1];
ugrid = g[[All, 1]]; weights = g[[All, 2]]; grid = Join[{0}, ugrid];

fd = NDSolveFiniteDifferenceDerivative[Derivative[1], grid]; m =
fd["DifferentiationMatrix"]; vart =
Table[u[i][t], {i, Length[grid]}]; varti =
Table[u[i]'[t], {i, 2, Length[grid]}]; vart1 =
Table[u[i]'[t], {i, Length[grid]}]; ux = m.vart; var =
Table[u[i], {i, Length[grid]}];

\[Mu][a_] := 1/(1 - a);
p0[x_] := If[x <= 1/2, (1 - 2 x)^3 (1 - x), 31 (2 x - 1)^3 (1 - x)];
u0[a_] := p0[a];

f[a_] := 2; int = Table[f[ugrid[[i]]] weights[[i]], {i, np}]; eqns =
Join[{D[u[1][t], t] == int.varti},
Flatten[Table[
u[i]'[t] == -ux[[i]] - \[Mu][grid[[i]]] u[i][t], {i, 2,
Length[grid]}]]];
ics = Table[u[i][0] == u0[grid[[i]]], {i, Length[grid]}];
tmax = 1;
{bb, mm} = CoefficientArrays[eqns, vart1];
rhs = -Inverse[mm].bb;
sol1 = NDSolve[{Table[vart1[[i]] == rhs[[i]], {i, Length[vart1]}],
ics}, var, {t, 0, tmax}];


With this code we calculate absolute error in Example 1 for np=40(left picture) and np=80 (right picture). For this code error is of $$h^{5/2}$$.

• I'm not sure you can take the analytical solution derived without the integral boundary condition, and use it to determine the boundary condition in NDSolve. You can calculate the basic reproduction number $R_0$ (expected number of offspring produced over a lifetime) for this model as NIntegrate[E^(-1/80*a) f[a], {a, 0, 100}] which is 1.672, so the population should grow, not go extinct. Commented Sep 20, 2020 at 2:14
• @ChrisK I used general solution of PDE with initial data $u(a,0)=u_0(a)$. This kind of solution exist and DSolve[] found it. The boundary condition $u(0,t)$ been used for numerical solution restricted on a. Commented Sep 20, 2020 at 11:12
• @ChrisK See update to my answer with comparison of your and my solutions. Commented Sep 20, 2020 at 13:05
• This isn't OP's original problem though. You decreased the birth rate f[a] from 0.1 to 0.025, so it's not surprising that the population dies off. The boundary condition represents reproduction, so it's an essential part of this model that can't just be changed to avoid NDSolve warnings. Commented Sep 20, 2020 at 14:35
• @ChrisK But this is unique solution with boundary and initial conditions are consistent. Also there are some not unique solutions with inconsistent boundary and initial conditions and you shown some of them. Commented Sep 20, 2020 at 15:40

To give some convincing publication to the round of answerers and the owner of the question: Numerical methods for the Lotka–McKendrick’s equation Galena Pelovska, Mimmo Iannelli∗ Dipartimento di Matematica, Universita degli Studi di Trento, via Sommarive 14, I-38050 Povo(Trento), Italy.

• Thanks, looks interesting Commented Sep 27, 2020 at 21:19