# How to solve $\frac{\partial}{\partial t} f(i,t) =b h(i) \int_0^\infty \mathrm {d} i\, h(i) f(i,t)$ with Mathematica?

Given $$b$$, $$h(i)$$, how do I solve the following in Mathematica for initial value $$f(i,0)=1$$? $$$$\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}$$$$

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

ClearAll["Global*"];
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}]


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

ClearAll["Global*"];
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}]


## 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.

ClearAll["Global*"];
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}];
ListContourPlot[Re@sol]
ListPlot[{Re@Last@sol, h$}, MultiaxisArrangement -> All, Filling -> Bottom]  • 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$. – josh Commented May 1, 2023 at 21:35 • @josh actually that appears to be the same, since 3 is$\frac{1}{h_2}$. How did you solve it? Commented May 2, 2023 at 0:00 • 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)$Commented May 2, 2023 at 0:12 • Ok thanks. Didn't see that. I posted how I did it below. – josh Commented May 2, 2023 at 12:18 ## 3 Answers 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).$$ Verification: $$\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}.$$ • 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 Commented May 1, 2023 at 23:59 • @YaroslavBulatov Yes, you are right, the constant$c$can be eliminated from the final result, see my edit. Commented May 2, 2023 at 8:44 • neat! I wonder if this method can also be used to solve mathematica.stackexchange.com/questions/284257/… Commented May 2, 2023 at 10:38 • @YaroslavBulatov I tried, but so far no progress... Commented May 2, 2023 at 11:11 • 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 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'} ]  with 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 == 1/2!Integral[f[x],{x}]Integral[f[y],{y}] =1/2!Integral[f[x]f[y],{x,y}]]  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}]  by f[x] == Integrate[Exp[I k x]g[k],{k}]  • 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)\$? 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}]


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