I would like to solve numerically the following integro-differential equation $$ \partial_t \rho(t,x) \,=\, \partial_x\big(f'(x)\,\rho(t,x)\big) \int_0^\infty f(\xi)\,\rho(t,\xi)\,d\xi \;+\\ +\; \partial_x\big(g'(x)\,\rho(t,x)\big) \int_0^\infty g(\xi)\,\rho(t,\xi)\,d\xi $$ where:
- $\rho$ is a probability distribution on $[0,\infty)$ which actually can degenerate to a convex combination of a Dirac delta and a density function;
- the initial condition $\rho(0,x)$ can be suitably chosen, such that $\int_0^\infty\rho(0,x)\,dx=1$;
- let's say the functions $f,g$ are given. They are strictly increasing, smooth but not analytic at $0$ indeed $f^{(k)}(0)=g^{(k)}(0)=0$ for all $k\geq1$.
I've tried with DSolve, but an exact solution is not found. Then I've tried with NDSolve and I get the following error:
NDSolve::delpde: Delay partial differential equations are not currently supported by NDSolve.
Is it possible to solve this equation using Mathematica? I am using Mathematica 11.
Edit
Here's the definition of $f,g$. Let $L(x)$ be a piecewise linear function taking value $l_0$ for $x\leq x_0$, $l_0+\frac{x-x_0}{x_1-x_0}\,(x_1-x_0)$ for $x_0\leq x\leq x_1$ and $l_1$ for $x\geq x_1$. Then set: $$ E(x) = \int_{-\infty}^{\infty} L(xz)\, \frac{e^{-\frac{z^2}{2}}}{\sqrt{2\pi}}\, dz $$ finally fix $c$ positive, $\epsilon\in(0,1)$ and let $$ f(x) = c\,E\big((1+\epsilon)\,x\big)-c \quad,\quad g(x) = c\,E\big((1-\epsilon)\,x\big)+c\;. $$ E.g. fix $l_0=-2.5,\,l_1=7.5,\,x_0=0.5,\,x_1=1.5$ and $c=1,\,\epsilon=0.6\,$.
Edit 2
I've obtained a plot of the solution implementing the Numerical Method of Lines suggested by @bbgodfrey , but there are same issues for $x$ close to $0$. Here's the resulting plot, from two points of view:
Solution $\rho(t,r)$ obtained by the numerical method of lines. View 1
Solution $\rho(t,r)$ obtained by the numerical method of lines. View 2
It seems something happens around $t\approx0.5$. What are those streight lines? Is there a way to see clearly the appearance of a Delta function and distinguish it from numerical problems?
Here's my code:
n = 1000; rmax = 5; T = 2;
X = Table[rmax/n*(i - 1), {i, 1, n + 1}];
Rho[t_] := Table[Subscript[ρ, i][t], {i, 1, n + 1}];
F = Table[f[X[[i]] $MachinePrecision], {i, 1, n + 1}];
G = Table[g[X[[i]] $MachinePrecision], {i, 1, n + 1}];
DF = Table[Df[X[[i]] $MachinePrecision], {i, 1, n + 1}];
DG = Table[Dg[X[[i]] $MachinePrecision], {i, 1, n + 1}];
(* Initial condition *)
gamma[r_] := 1/(Gamma[k] θ^k) r^(k - 1) Exp[-r/θ]
k = 10; θ = 0.1;
ic = Thread[ Drop[Rho[0], -1] == Table[gamma[X[[i]]], {i, 1, n}] ];
(* Boundary condition *)
Subscript[ρ, n + 1][t_] := 0
(* ODE's *)
rhs[t_] :=
ListCorrelate[{-1, 1}, DF*Rho[t]]*Total[F*Rho[t]] +
ListCorrelate[{-1, 1}, DG*Rho[t]]*Total[G*Rho[t]]
lhs[t_] := Drop[D[Rho[t], t] , -1]
eqns = Thread[lhs[t] == rhs[t]];
lines =
NDSolve[
{eqns, ic}, Drop[Rho[t], -1], {t, 0, T},
Method -> {"EquationSimplification" -> "Residual"}];
ParametricPlot3D[
Evaluate[Table[{rmax/n*i, t, First[Subscript[ρ, i][t] /. lines]}, {i, 1, n/2}]],
{t, 0, 1},
AxesLabel -> {"r", "t", "ρ"}, BoxRatios -> {1, 1, 1}]
x
and provide boundary conditions. Even then, `NDSolve cannot solve this problem without help. Please see my solution to a related problem here. $\endgroup$x
to obtain an expression for ρ(t,0). Note, though, that you will encounter resolution problems as your solution collapses to a delta-function. Perhaps, representing rho as the sum of a numerical function to be calculated plus a one-parameter function peaking atx = 0
, with that parameter also to be calculated. Adding the definitions off[x]
andg[x]
might be helpful to readers. $\endgroup$