Using DSolve for integro-differential equation $\frac{\partial}{\partial t} f =hf+ \int \mathrm {d} i\, h f$?

Edit Apr 28 Wolfram Support confirmed that DSolve behavior is unexpected and escalated to developers, Case 5025062. Meanwhile, I'm interested in a work-around or an alternative way to solve this with Mathematica.

DSolve fails to give expected solution for integro-differential equation of the following form. $$\frac{\partial}{\partial t} f(t, i) =a h(i) f(t, i) + b h(i) \int \mathrm {d} i_1\, h(i_1) f(t, i_1)$$

IE, solving for $$f$$ using DSolve and then plugging the solution back into the equation gives False as seen below. Any suggestion how to find a working solution?

ClearAll["Global*"];
h = (i + 1)^-2;
inti[expr_] := Integrate[expr, {i, 0, Infinity}];
eqn = D[f[i, t], t] == h f[i, t] + h inti[h f[i, t]];
sol = DSolveValue[{eqn, f[i, 0] == 1}, f[i, t], {i, t}]
f0 = Function @@ {{i, t}, sol};

Print["Solution: ", f0[i, t]];
Print["Back-substitute: ", eqn0 = eqn /. f -> f0];
Print["Solution correct: ", eqn0 /. {i -> 2, t -> 3}]  (*False*)


This is the continuous analogue of the DPR1 problem described here.

NDSolve also doesn't give expected solution. It's solution appears identical to the one where you make substitution $$\int_0^\infty \mathrm{d}i\ h(i) f(t, i)=f(t, i)\int_0^\infty \mathrm{d}i\ h(i)$$. However, this transformation gives a different problem

ClearAll["Global*"];
h = (i + 1)^-2;
inti[expr_] := Integrate[expr, {i, 0, Infinity}];
eqn = D[f[i, t], t] == -2 h f[i, t] + h inti[h f[i, t]];
solSto =
NDSolveValue[{eqn, f[i, 0] == 1}, f[i, t], {t, 0, 10},
i \[Element] ImplicitRegion[i > 0 && i < 10, i]];

eqn = D[f[i, t], t] == -2 h f[i, t] + h f[i, t];
solDet =
NDSolveValue[{eqn, f[i, 0] == 1}, f[i, t], {t, 0, 10},
i \[Element] ImplicitRegion[i > 0 && i < 10, i]];

solDet == solSto (* true, but this shouldn't be *)
ContourPlot[solSto, {i, 0, 10}, {t, 0, 10}]


I can solve this differential equation numerically by discretizing, and solving discrete linear system using inverse Laplace transform, but that doesn't really scale past 10 dimensions

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 = -2 DiagonalMatrix[h$$] + {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[sol]


Below you can see solution of NDSolve vs the expected one (from Laplace Transform)

• I tried Block[{h = Exp[-i]}, ResourceFunction["TraceView"][solve[h], TraceInternal -> True]]; to navigate through the internal implementation but I did not see anything particularly instructive. Maybe I did not look well enough Apr 26, 2023 at 1:01
• An advantage of mathematical notations is that they allow to convey the statement unambiguously. But your writing is actually ambiguous. It is not clear what equation are you actually trying to solve. Why not write definite integral as definite integral, and use different symbols for variables in different context? There is a good proof for what I said, as MA itself misunderstood you. Apr 27, 2023 at 16:48
• @yarchik I'm not much of a mathematician, but the Mathematica code should make the problem unambiguous, no? I believe it's a bug in Mathematica, the solution I get from DSolve doesn't satisfy the differential equation given to DSolve, edited example of clarity Apr 27, 2023 at 17:55
• Could you expand on your post a bit better so that we can have some of the insights that your work has provided you so far? that might help us futher.
– alex
Apr 27, 2023 at 17:59
• As it stands, its is abuse of notation. The integration variable must not be used outside the integral. If the algebraical use int f(i) di = F(i) with F'=f is meant, extensive use of Derivative should be use, to clarify things: Derivative[0,1][f] = f + Derivative[-1,0] [ f[#1,#2] g[#1]&]. Translated to Integrate, include the limits and localize the summation variable i, di ->j , dj . To improve human readabiliy the dummy variables should resemble the global names. Apr 28, 2023 at 23:11

modified Infinite integration range {z,0,Infinity} has to be numerically adapted!

case h=1/(1+x)^2,a=-2,b=1

h=1/(1+#)^2&;
a=-2;
b=1;


Numerical solution(NDSolve, iteration NestList or FixedPoint)

inf = 10000; (*~Infinity *)
nl = NestList[Function[{fa}, Block[{int, x, t },
int[t_?NumericQ]:=Block[{z}, NIntegrate[ h[z] fa[z, t], {z, 0, inf(* Infinity*) },Method -> {Automatic, "SymbolicProcessing" -> 0}]];
NDSolveValue[{Derivative[0, 1][f][x,t] ==a h[x] f[x, t] + b h[x] int[t]   , f[x, 0] == 1},
f, {x, 0, inf }, {t, 0, 10}, MaxStepSize -> {100, Automatic} ] ]],
1 &(*Identity*), 5] ;
Plot3D[nl[[-1]][x, t], {x, 0, 100}, {t, 0, 10}, PlotRange -> All,
MeshFunctions -> {#3 &}, AxesLabel -> {x, t, f[x, t]}]


final modification (14.5.2023)

The iterative numerical solution might be improved considerably with the tranformation u->h[x], f[x,t]->fu[u=h[x],t]

The integrodifferential equation transforms to

Numerical iteration now converges much faster.

nlu = NestList[
Function[{fa},
Block[{int, x, t},int[t_?NumericQ] :=
Block[{z}, NIntegrate[ 1/(2 Sqrt[z]) fa[z, t], {z, 0, 1 },  Method -> {Automatic, "SymbolicProcessing" -> 0},PrecisionGoal -> 5,AccuracyGoal -> 10]];
NDSolveValue[{Derivative[0, 1][f][u, t] == a u f[u, t] + b u int[t], f[u, 0] == 1}, f, {u, 0, 1}, {t, 0, 10} ]]], 1 &(*Identity*), 10]; // AbsoluteTiming (* 1.9s*)

Plot3D[nlu[[-1]][h[x] , t], {x, 0, 10 }, {t, 0, 10}, PlotRange ->All,MeshFunctions -> {#3 &}, AxesLabel -> {x, t, f[x, t]},PlotPoints-> 100]MeshFunctions -> {#3 &}, AxesLabel -> {x, t, f[x, t]},PlotPoints -> 100]


This solution is nearly identical to the solution found by Laplacetransformation !

end of final modification (14.5.2023)

Analytical solution(powerseries)

Integrating the integrodifferentialequation vs time , thereby considering the ic f[x,0]==1, gives

fixedpoint equation

F = Function[{x, t},
Evaluate[1 + a h[x] Integrate[  #[x, \[Tau]] , {\[Tau], 0, t}] + b h[x]Integrate[h[z] #[z, \[Tau]], {z, 0, Infinity}, {\[Tau], 0, t}]]] &

fp = FixedPoint[F, 1 &, 5]
Collect[fp[x, t], t, Expand]
(*
1 - t/(1 + x)^2 + t^2 (1/(1 + x)^4 - 1/(6 (1 + x)^2)) +
t^4 (1/(3 (1 + x)^8) - 1/(18 (1 + x)^6) - 13/(540 (1 + x)^4) - 323/(
22680 (1 + x)^2)) +
t^5 (-(2/(15 (1 + x)^10)) + 1/(45 (1 + x)^8) + 13/(1350 (1 + x)^6) +
323/(56700 (1 + x)^4) + 6647/(1701000 (1 + x)^2)) +
t^3 (-(2/(3 (1 + x)^6)) + 1/(9 (1 + x)^4) + 13/(270 (1 + x)^2))
*)


From this result one might guess the general form of the solution f[x,t]

The coefficients c[k,j] follow from equation null==0 comparing powers of t

solF = Function[{x, t}, 1 + Sum[t^k c[k, j] (1 + x)^(-2  j ), {k, 1,n}, {j, 1, k }]]
null = Collect[1 + a h[x] Integrate[  solF[x, \[Tau]] , {\[Tau], 0, t}] +b h[x] Integrate[h[z] solF[z, \[Tau]], {z, 0, Infinity}, {\[Tau], 0, t}] - solF[x, t], t, FullSimplify]
(*-((t (1 + c[1, 1]))/(1 + x)^2) +
(t^2 ((-5 + x (2 + x)) c[1, 1] -6 ((1 + x)^2 c[2, 1] + c[2, 2])))/(6 (1 + x)^4) +
...*)


{c[1, 1] -> -1, c[2, 1] -> -(1/6), c[2, 2] -> 1,...} 

• Very interesting! (btw, the interesting cases are a=-2, b=1 and a=-2, b=0, they correspond to fully stochastic vs. fully deterministic version of DPR1 problem May 3, 2023 at 8:56
• Ok I' ll modify my answer to the case a=-2,b=1! May 3, 2023 at 8:58
• thanks, will dig into it. It's interesting that you integrate over time rather than take Laplace Transform. For the discrete version of the problem I can apply Sherman-Woodbury to get a simple formula for laplace transform of $g(t)=\int \mathrm{d}i\ f(i, t)$, $$\mathcal{L}(g)(s)=\langle 1/z,1\rangle + \frac{b}{1-b\langle h/z,h\rangle} \langle h/z,1\rangle^2\\ z=s-a h$$ May 3, 2023 at 9:25
• In your first solf, why do you call NDSolve as a subroutine to fixed point iteration rather than just using NDSolve directly? May 3, 2023 at 10:07
• It's an iterative solution where I evaluate the integralpart of the ode with the last iteration. The solutions nl should converge against the unknown solution(hopefully) May 3, 2023 at 10:21

Is it possible, you try to solve someting like the damped linear equation with

 f'[t] == - λ  f[t] + h[t] f[t]  ->  Log[f[t]]] == \int^t( h[s] -λ ds

f[t] == f[0] Exp[-λ t +\int_0^t h[s]ds ]


Under integral transforms the product maps to a convolution integral. I had the suspect h as 1/(I+s) that is the Laplace transform of e^-t

• As far as I can say, f` is a function of two arguments. In your solution it is a function of just one. May 14, 2023 at 13:46