5
$\begingroup$

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}]
Improve this question
$\endgroup$
15
  • 1
    $\begingroup$ In general, a differential equation cannot be solved numerically over an infinite domain, so you will need to truncate the problem in x and provide boundary conditions. Even then, `NDSolve cannot solve this problem without help. Please see my solution to a related problem here. $\endgroup$ Commented Sep 18, 2020 at 20:26
  • $\begingroup$ @bbgodfrey thank you, I'm trying to understand this numerical method of lines. Before I'd like to ask if it's possible to avoid boundary conditions. Actually in my problem I have the constraint $\int_0^\infty\rho(t,x)dx=1$ for all $t\geq0$. I could impose $\rho(t,\bar x)=0$ where $\bar x=100$ (or any suitably large number), but I have no condition to impose for $\rho(t,0)$. Furthermore from other results I expect $\rho(t,0)$ to change a lot during time passing from $0$ to a Dirac delta $\endgroup$ Commented Sep 20, 2020 at 14:57
  • $\begingroup$ It is not possible, even in principle, to avoid boundary conditions. However, you could try integrating the entire integral-differential equation over 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 at x = 0, with that parameter also to be calculated. Adding the definitions of f[x] and g[x] might be helpful to readers. $\endgroup$ Commented Sep 20, 2020 at 16:32
  • 1
    $\begingroup$ So much for that idea of mine. Adding the definitions of f[x] and g[x] to the question would be helpful. $\endgroup$ Commented Sep 20, 2020 at 18:21
  • 1
    $\begingroup$ In your 1st block of code, when making tables, you scale various functions by $MachinePrecision. I've find that a strange choice. Can you tell why you chose that scale factor? $\endgroup$ Commented Sep 24, 2020 at 23:39

2 Answers 2

Reset to default
4
$\begingroup$

Since in the original code there are instabilities due to low order approximation we can use 4th order numerical algorithm I have developed for the Lotka-McKendrick demographic model (see the very last code in my answer). First we define function f, g using next exact expression for $E(x)$:

l0 = -25/10; l1 = 75/10; x0 = 1/2; x1 = 3/2; c = 1; eps = 3/5; L[x_] := Piecewise[{{l0, x <= x0}, {l0 + (l1 - l0) (x - x0)/(x1 - x0), x0 < x <= x1}, {l1, x > x1}}]; Integrate[L[x z] Exp[-z^2/2], {z, -Infinity, Infinity}, Assumptions -> {x > 0}]/Sqrt[2 Pi] (*1/(4 Sqrt[2 \[Pi]])5 \[ExponentialE]^(-(9/(8 x^2))) (-\ \[ExponentialE]^((9/(8 x^2))) Sqrt[2 \[Pi]]-8 x+8 \ \[ExponentialE]^(1/x^2) x+2 \[ExponentialE]^(9/(8 x^2)) Sqrt[2 \[Pi]] \ Erf[1/(2 Sqrt[2] x)]-3 \[ExponentialE]^(9/(8 x^2)) Sqrt[2 \[Pi]] \ Erf[3/(2 Sqrt[2] x)]+3 \[ExponentialE]^(9/(8 x^2)) Sqrt[2 \[Pi]] \ Erfc[3/(2 Sqrt[2] x)])*)

Therefore we can explicit define functions $f(x),g(x),E(x),E'(x)f'(x), g'(x)$ as f,g,eL,eL1,df,dg, we have

eL[x_] := 1/(4 Sqrt[2 \[Pi]]) 5 E^(-(9/( 8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]); eL1[x_] := ( 45 E^(-(9/( 8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]))/( 16 Sqrt[2 \[Pi]] x^3) + ( 5 E^(-(9/( 8 x^2))) (-8 + 8 E^(1/x^2) + (9 E^(9/(8 x^2)) Sqrt[\[Pi]/2])/( 2 x^3) + 18/x^2 - (18 E^(1/x^2))/x^2 - ( 9 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[1/(2 Sqrt[2] x)])/x^3 + ( 27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[3/(2 Sqrt[2] x)])/(2 x^3) - ( 27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erfc[3/(2 Sqrt[2] x)])/(2 x^3)))/( 4 Sqrt[2 \[Pi]]); f[x_] := c eL[(1 + eps) x] - c; df[x_] := c (1 + eps) eL1[(1 + eps) x]; g[x_] := c eL[(1 - eps) x] + c; dg[x_] := c (1 - eps) eL1[(1 - eps) x];

Second step, we call

Needs["DifferentialEquations`NDSolveProblems`"]; Needs["DifferentialEquations`NDSolveUtilities`"]; Get["NumericalDifferentialEquationAnalysis`"];

Now we define grid and weights for numerical integration using GaussianQuadratureWeights[] and DifferentiationMatrix on the same grid using FiniteDifferenceDerivative:

np = 100; gqw = GaussianQuadratureWeights[np, 0, 5]; ugrid = gqw[[All, 1]]; weights = gqw[[All, 2]]; fd = NDSolve`FiniteDifferenceDerivative[Derivative[1], ugrid]; m = fd["DifferentiationMatrix"];

Finally we define all needed vectors, matrixes, equations and solve system of ODEs using NDSolve

Quiet[varf = Table[df[ugrid[[i]]] u[i][t], {i, Length[ugrid]}]; varg = Table[dg[ugrid[[i]]] u[i][t], {i, Length[ugrid]}]; varu = Table[u[i][t], {i, Length[ugrid]}]; var = Table[u[i], {i, Length[ugrid]}]; ufx = m.varf; ugx = m.varg; intf = Table[f[ugrid[[i]]] weights[[i]], {i, np}]; intg = Table[g[ugrid[[i]]] weights[[i]], {i, np}]]; u0[r_] := 1/(Gamma[k] \[Theta]^k) r^(k - 1) Exp[-r/\[Theta]] k = 10; \[Theta] = 0.1; ics = Table[u[i][0] == u0[ugrid[[i]]], {i, np}]; eqns = Table[D[u[i][t], t] == ufx[[i]] (intf.varu) + ugx[[i]] (intg.varu), {i, np}]; tmax = 2; sol = NDSolve[{eqns, ics}, var, {t, 0, tmax}, Method -> {"EquationSimplification" -> "Residual"}];

Visualization of numerical solution

lst = Flatten[ Table[{t, ugrid[[i]], u[i][t] /. sol[[1]]}, {t, 0, 2, 1/50}, {i, np}], 1]; ListPlot3D[lst, Mesh -> None, PlotRange -> All, AxesLabel -> {"t", "x"}]

Figure 1 We can compare this result with the original code running for n=50(left picture) and n=100 (right). On the left picture we can recognize solution shown above. But there are also unphysical oscillation with amplitude increasing 10 times with n increases from 50 to 100. Original code as I am used for n=50

eL[x_] := 1/(4 Sqrt[2 \[Pi]]) 5 E^(-(9/( 8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]); eL1[x_] := ( 45 E^(-(9/( 8 x^2))) (-E^((9/(8 x^2))) Sqrt[2 \[Pi]] - 8 x + 8 E^(1/x^2) x + 2 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[1/(2 Sqrt[2] x)] - 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erf[3/(2 Sqrt[2] x)] + 3 E^(9/(8 x^2)) Sqrt[2 \[Pi]] Erfc[3/(2 Sqrt[2] x)]))/( 16 Sqrt[2 \[Pi]] x^3) + ( 5 E^(-(9/( 8 x^2))) (-8 + 8 E^(1/x^2) + (9 E^(9/(8 x^2)) Sqrt[\[Pi]/2])/( 2 x^3) + 18/x^2 - (18 E^(1/x^2))/x^2 - ( 9 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[1/(2 Sqrt[2] x)])/x^3 + ( 27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erf[3/(2 Sqrt[2] x)])/(2 x^3) - ( 27 E^(9/(8 x^2)) Sqrt[\[Pi]/2] Erfc[3/(2 Sqrt[2] x)])/(2 x^3)))/( 4 Sqrt[2 \[Pi]]); f[x_] := c eL[(1 + eps) x] - c; df[x_] := c (1 + eps) eL1[(1 + eps) x]; g[x_] := c eL[(1 - eps) x] + c; dg[x_] := c (1 - eps) eL1[(1 - eps) x]; n = 50; rmax = 5; T = 2; X = Table[rmax/n*(i - 1) + 10^-6, {i, 1, n + 1}]; Rho[t_] := Table[Subscript[\[Rho], i][t], {i, 1, n + 1}]; F = Table[f[X[[i]] ], {i, 1, n + 1}]; G = Table[g[X[[i]] ], {i, 1, n + 1}]; DF = Table[df[X[[i]]], {i, 1, n + 1}]; DG = Table[dg[X[[i]] ], {i, 1, n + 1}]; (*Initial condition*) gamma[r_] := 1/(Gamma[k] \[Theta]^k) r^(k - 1) Exp[-r/\[Theta]] k = 10; \[Theta] = 0.1; ic = Thread[Drop[Rho[0], -1] == Table[gamma[X[[i]]], {i, 1, n}]]; (*Boundary condition*) Subscript[\[Rho], 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"}];

Visualization of numerical solutions for n=50 (left) and n=100 (right)

lst = Table[{t, X[[i]], Subscript[\[Rho], i][t] /. lines[[1]]}, {t, 0, T, 1/25}, {i, n}]; ListPlot3D[Flatten[lst, 1], ColorFunction -> "Rainbow", Mesh -> None, AxesLabel -> {"t", "x", ""}, PlotRange -> All]

Figure 2

Improve this answer
$\endgroup$
1
$\begingroup$

This is not an answer, but some comments on solving this type of problem that are too long and to made in comments to the question.

Regarding scaling up and down: In my opinion, in order to become proficient in solving difficult problems it is imperative to learn how to scale the problem down and then back up again. For example, you have: $$ \frac{\partial \rho}{\partial t}=\frac{\partial}{\partial t}\left(f'\rho\right)\int_0^{\infty} f(x)\rho(t,x)dx+\cdots $$ Notice the dots. When removed, that's scalling it down to a simpler form. Can you solve just that one? Maybe though it has no solution. I don't know. How about take out the $f'\rho$ term, say:

$$ \frac{\partial \rho}{\partial t}+\frac{\partial p}{\partial x}=\int_0^{\infty} f(x)\rho(t,x)dx $$

That one? How about take out the $f(x)$ term in the integrand? Just how much do you have to scale it down while still retaining it's PIDE nature in order to solve it? How about just solving any simple (similar somewhat) PIDE to perfect the method and then add complexity (terms) to the problem until you reach the equation you wish to solve.

Of course that takes a lot of work and sometimes you will of course run into problems where scaling it up further causes a significant hurtle to solve. But surprisingly, this method has frequently been very successful with tough problems I have worked on but not always. Here's an example: $$ f+\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}=\int_x^{\infty} \int_y^{\infty}f(u,v)dudv $$ beautiful huh, but a bit intimidating. How about we scale it down: $$ f+\frac{df}{dx}=\int_x^{\infty} f(u)du $$ That's easier and it turns out, the solution to that one leads easily to the solution of the first. :)

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for your suggestion, but I don't think this is an answer (by the way I didn't downvote) $\endgroup$ Commented Sep 24, 2020 at 11:35

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.