Given $b$, $h(i)$, how do I solve the following in Mathematica for initial value $f(i,0)=1$? $$\begin{equation}\frac{\partial}{\partial t} f(i,t) =b h(i) \int_0^\infty \mathrm {d} i\, h(i) f(i,t)\label{eq}\tag{1}\end{equation}$$

An internet user says that solution is $f(i,t)=1+\frac{1}{h_2} (e^{bh_2 t} - 1)h(i)$. I suspect it is incorrect.

Any tips?

What I tried:

DSolve: returns unevaluated

Let $h(i)=(1+i)^{-2}, b=1$, we can use DSolve

h = (i + 1)^-2;
inti[expr_] := Integrate[expr, {i, 0, Infinity}];
eqn = D[f[i, t], t] == h inti[h f[i, t]];
DSolveValue[{eqn, f[i, 0] == 1}, f[i, t], {t, i}]

enter image description here

Meanwhile, the following DSolve succeeds but the output doesn't satisfy \eqref{eq}, similar to issue here

DSolveValue[{eqn, f[i, 0] == 1}, f[i, t], {i, t}]

NDSolve: produces numerical noise

h = (i + 1)^-2;
inti[expr_] := Integrate[expr, {i, 0, Infinity}];
eqn = D[f[i, t], t] == h inti[h f[i, t]];

numSteps = 100;
maxDim = 10;
solSto = 
  NDSolveValue[{eqn, f[i, 0] == 1}, f[i, t], {t, 0, numSteps}, 
   i \[Element] ImplicitRegion[i >= 0 && i <= maxDim, i]];
ContourPlot[solSto, {i, 0, maxDim}, {t, 0, numSteps}]

func[i0_] := Block[{i = i0, t = numSteps}, solSto]
Plot[{h, func[i]}, {i, 0, 9}]

enter image description here

InverseLaplaceTransform: slow

Discretizing at $n$ points and applying inverse Laplace transform gives something that's within a constant factor to the solution provided. However it

  1. Doesn't give me the missing constant factor

  2. Is very slow. I believe it does $O(n^3)$ matrix inversion, and then issues $O(n^2)$ calls to inverse Laplace transform.

d = 10;
h = (i + 1)^-2;
discretize[expr_] := Table[N@expr /. {i -> i0}, {i0, 0, d - 1}];
{f0$, h$} = discretize /@ {1, h};

A = {h$}\[Transpose] . {h$};
ii = IdentityMatrix[d];

f$ = Function @@ {{t}, 
    InverseLaplaceTransform[Inverse[s*ii - A], s, t] . f0$};
sol = Table[f$[t], {t, 0, 10}];
ListPlot[{Re@Last@sol, h$}, MultiaxisArrangement -> All, 
 Filling -> Bottom]

enter image description here

  • $\begingroup$ Using your reference in the Mathematica forum and letting $f(i,t)=1+g(t)h(i)$ I get a different answer than above: $$\displaystyle f(i,t)=1+\frac{1}{(1+i)^2}\left(3 e^{b h_2 t}-\frac{1}{h_2}\right)$$ with $\displaystyle h_2=\int_0^{\infty} h(i)^2 di=1/3$ and this solution satisfies the IDE and the initial conditions $f(i,0)=1$. $\endgroup$
    – josh
    Commented May 1, 2023 at 21:35
  • $\begingroup$ @josh actually that appears to be the same, since 3 is $\frac{1}{h_2}$. How did you solve it? $\endgroup$ Commented May 2, 2023 at 0:00
  • $\begingroup$ The solution I get by heuristically extrapolating from discrete time version of the problem is $$f(i, t)=1+\frac{h_1}{h_2}(e^{bh_2 t} - 1)$$ where $h_1=\int_0^\infty \mathrm{d}i\ h(i)$ $\endgroup$ Commented May 2, 2023 at 0:12
  • $\begingroup$ Ok thanks. Didn't see that. I posted how I did it below. $\endgroup$
    – josh
    Commented May 2, 2023 at 12:18

3 Answers 3


One possible solution of

$$\frac{\partial}{\partial t} f(x,t) =b h(x) \int_0^\infty \mathrm {d} z\, h(z) f(z,t)$$

could be obtained using the Ansatz $f(x,t)=g(x)+h(x)F(t)$. We have $$ \frac{d}{d t}F(t)=b \int_0^\infty \mathrm {d} z\, h(z)^2F(t)+b\int_0^\infty \mathrm {d} z\, h(z)g(z). $$ Denote $$ H=b \int_0^\infty \mathrm {d} z\, h(z)^2,\\ G=b\int_0^\infty \mathrm {d} z\, h(z)g(z). $$ and solve the resulting ODE for F(t) $$ \frac{d}{d t}F(t)=HF(t)+G. $$ What are the boundary conditions? From $𝑓(x,0)=1$ we infer $$ g(x)+h(x)F(0)=1. $$ From this equation follows $$ F(0)=\frac{1-g(x)}{h(x)}. $$ Thus, $g(x)=1-ch(x)$ and $F(0)=c$. Now

DSolveValue[{D[F[t], t] == H F[t] + G, F[0] == c}, F[t], t]

leading to $$ F(t)=\frac{c H e^{H t}+G e^{H t}-G}{H}. $$ Also, per definition $$ G=b\int_0^\infty \mathrm {d} z\, h(z)(1-ch(z))=K-cH,\\ K=b\int_0^\infty \mathrm {d} z\, h(z). $$

The full solution is then $$ f(x,t)=1-\frac{K}{H}h(x)+h(x)\frac{K}{H} e^{H t},\\ H=b \int_0^\infty \mathrm {d} z\, h(z)^2,\\ K=b\int_0^\infty \mathrm {d} z\, h(z). $$


$$ \text{lhs}=\frac{\partial}{\partial t} f(x,t)=h(x)K e^{H t}.\\ b \int_0^\infty \mathrm {d} z\, h(z) f(z,t)= b \int_0^\infty \mathrm {d} z\, h(z) \left(1-\frac{K}{H}h(z)+h(z)\frac{K}{H} e^{H t}\right)=K-\frac{K}{H}H+\frac{HK}{H}e^{H t}=Ke^{H t},\\ \text{rhs}=b h(x)\int_0^\infty \mathrm {d} z\, h(z) f(z,t)=h(x)Ke^{H t},\\ \text{lhs}=\text{rhs}. $$

  • $\begingroup$ It's weird that there's not enough boundary conditions given that the discrete version of this problem is $f_{t+1}-f_t=bhh^T f_t$ which is uniquely defined by value of $f_0$. We can define corresponding problem for finer discretization, which will also be unique, hence I expect the continuous limit to be unique $\endgroup$ Commented May 1, 2023 at 23:59
  • $\begingroup$ @YaroslavBulatov Yes, you are right, the constant $c$ can be eliminated from the final result, see my edit. $\endgroup$
    – yarchik
    Commented May 2, 2023 at 8:44
  • $\begingroup$ neat! I wonder if this method can also be used to solve mathematica.stackexchange.com/questions/284257/… $\endgroup$ Commented May 2, 2023 at 10:38
  • $\begingroup$ @YaroslavBulatov I tried, but so far no progress... $\endgroup$
    – yarchik
    Commented May 2, 2023 at 11:11
  • $\begingroup$ I made some progress solving the discrete version for $h(t)=\int_0^\infty \mathrm{d}i\ f(i,t)$. Using Sherman-Morrison formula, I verified the following: $$\mathcal{L}(g)(s)=\langle 1/z,1\rangle + \frac{b}{1-b\langle h/z,h\rangle} \langle h/z,1\rangle^2$$ where $z=s-ah$, $1$ refers to vector of ones, $x/y$ is pointwise division of vectors $\endgroup$ Commented May 2, 2023 at 11:53

Seems to be a the standard problem of a parabolic evolution equations, depending on the integral kernel h.

Your discrete comment's definition translates to

 f[t+1,n]-f[t,n] = b Sum[ Sum[ h[n,m] h[k,m], {m} ]  f[t,k],{k}] 

(don't be sparing at spaces, at least not at implicit "*"'s)

The continuous t,k limit is

dt d/dt f[t,k] = b Integrate[#,{k'}]& [Integrate[ h[k',k''] h[k'',k'] f[t,k'], {k''}]]  
== Integrate[  K[k,k'] f[t,k'], {k'}  ]


K[k,k'] == K[k',k] 

a symmetric kernel.

A symmetric kernel K defines a selfadjoint linear operator A, written as a linear function

f' =  A [f]  ->  f[t] = Exp[A t][f[0]]

if A has negative spectrum only, the matrix Exp converges for all t > 0. It's of course a delicate problem because nth-powers of integrals are in fact n-dimensional integrals over powers of f

    1/2! Integral[f[x],{x}]^2 == 

Eg if

K = delta(k-k') k'^2  

 f'[t][k] = -k^2 f[t][k]

that is the Fourier transform of

g'[t][x]  = D[g[t][x],{x,2}]


f[x] == Integrate[Exp[I k x]g[k],{k}]
  • $\begingroup$ Interesting....I don't quite get the point of $g$ in the last 2 equations, is that the continuous analogue of formula $f(t)=\mathcal{L}^{-1}(s I - A)$ in the discrete-time evolution $f(t+1)-f(t)=A f(t)$? $\endgroup$ Commented May 2, 2023 at 7:45

Assuming solution $f(i,t)=1+g(t)h(i)$ and substituting into the IDE: $$ \begin{align} h(i)\frac{dg}{dt}&=b h(i)\int_0^{\infty} h(i)\left(1+g(t)h(i)\right)di \\ &=b h(i)\left(\int_0^{\infty} h(i)di+g(t)\int_0^{\infty} (h(i))^2di\right)\\ &=b h(i)\left(1+g(t)h_2\right) \end{align} $$ So that $$\displaystyle h(i)g'(t)=b h(i)\left(1+g(t)h_2\right)$$ Solving via Mathematica code:

theSol = DSolve[h[i] g'[t] == b h[i] (1 + g[t] h2), g[t], t][[1]]
gSol[t_] = g[t] /. theSol
(* -(1/h2) + E^(b h2 t) C[1] *)

then $\displaystyle g(t)=-\frac{1}{h2}+c_1 e^{bh_2 t}$ and therefore: $$ f(i,t)=1+h(i)\left(-\frac{1}{h_2}+c_1 e^{bh_2 t}\right) $$ and since $h_2=1/3$, we let $c_1=3$ since $f(i,0)=1$ or $$ f(i,t)=1+\frac{1}{(1+i)^2}\left(3 e^{bh_2 t}-3\right) $$ Plotting the solution for $b=2$:

b = 2;
h2 = 1/3;
fSol[i_, t_] := 1 + (-3 + 3 E^(b h2 t))/(1 + i)^2;
Plot3D[fSol[i, t], {i, 0, 5}, {t, 0, 5}, 
 PlotRange -> All, 
 AxesLabel -> {Style["i", 14, Bold, Black], 
   Style["t", 14, Bold, Black], Style["f", 14, Bold, Black]}, 
 BoxRatios -> {1, 1, 1}]

enter image description here

  • $\begingroup$ Thanks, it seems pretty clear NDSolve does not give correct result here, contacted support $\endgroup$ Commented May 2, 2023 at 13:21
  • $\begingroup$ Compared NDSolve to this solution here, looks like this. May be missing some secret flags to make NDSolve work $\endgroup$ Commented May 2, 2023 at 16:17

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.