# Numerically solving a multivariable integro-differential equation

I have an integro-differential equation of the form,

$$\small\frac{\partial f(x,t)}{\partial t}=\int_{-5}^5 |x-y|\,f\left(-\frac{x}{3}+\frac{4y}{3},t\right)\,f\left(\frac{2x}{3}+\frac{y}{3},t\right)\,\mathrm dy-\int_{-5}^5 |x-y|\,f(x,t)\,f(y,t)\,\mathrm dy$$

I have the function $$f(x,0)$$ representing the function at $$t=0$$. I want to numerically integrate the above equation to find the function $$f(x,t)$$ at any time $$t$$. I tried using NDSolveValue but Mathematica says the system is underdetermined. I'll be glad if someone can let me know if this can be solved using Mathematica.

This is my code, FYI

x1[x_, y_] = -x/3 + 4*(y/3);

y1[x_, y_] = 2*(x/3) + y/3;

eqn = D[f[x, t], t] == Integrate[Abs[x - y]*f[x1[x_, y_], t]*f[y1[x_, y_], t], {y, -5, 5}] - Integrate[Abs[x - y]*f[x, t]*f[y, t], {y, -5, 5}];

init = f[x, 0] == 1/10;

sol2 = NDSolveValue[{eqn, init}, f[x, t], {{x,-5,5}, {t, 0, 100}}];


Thanks a lot,

• From Documentation NDSolve or NDSolveValue this type equation can't solve. – Mariusz Iwaniuk Apr 12 at 16:34
• You equation does not contain any spatial derivatives and no spatial singularities, and so you may be able to solve it using a simple interpolation scheme. Even order-zero interpolation may be sufficient here if the grid is chosen fine enough. – Roman Apr 13 at 11:21
• The solution is the trivial $f(x,t)=\frac{1}{10}$, isn't it? – xzczd Apr 14 at 10:10
• No..f(x,t) behaves like a gaussian at large times...The above equation is a form of the famous Boltzmann transport equation, which converges to maxwell distribution at large times. – sammy Apr 14 at 11:08
• You need to add @xzczd in your comment or I won't get the reminder. Why do you think $f(x,t)=\frac{1}{10}$ isn't a solution? – xzczd Apr 14 at 11:20

We need a functional to evaluate the right-hand-side from a given function $$f$$:

xmin = -5;
xmax = 5;
rhs[f_, x_ /; xmin <= x <= xmax] :=
NIntegrate[Abs[x - y]*f[(4 y - x)/3]*f[(2 x + y)/3], {y, xmin, xmax}] -
f[x]*NIntegrate[Abs[x - y]*f[y], {y, xmin, xmax}]


Given a function $$f_t(x)$$ and a time step $$\Delta t$$, propagate the function forward in time (using the Euler method) and return the new function $$f_{t+\Delta t}(x)$$:

propagate[f_, Δt_?Positive] :=
Evaluate@Piecewise[{{FunctionInterpolation[f[x] + Δt*rhs[f, x], {x, xmin, xmax}][#],
xmin <= # <= xmax}}] &


The integro-differential equation can now be solved: starting from a function F[0] at $$t=0$$ and progressing by time-steps $$\Delta t$$ using the simple Euler method above, using partial memoization to save an InterpolatingFunction object at every time step:

Δt = 1/100;
Clear[F];
F[0] = Piecewise[{{1/(xmax - xmin), xmin <= # <= xmax}}] &;
F[t_?Positive] := F[t] = propagate[F[t - Δt], Δt]
F[x_?NumericQ, t_?NonNegative] := F[Round[t, Δt]][x]


The function F[x,t] represents the solution, but keep in mind that the $$t$$-coordinate has a resolution of $$\Delta t$$ and the $$x$$-coordinate has a finite resolution because it is represented as an interpolation. Using a smaller $$\Delta t$$ makes the solution more accurate but takes more time.

ContourPlot[F[x, t], {t, 0, 2}, {x, xmin, xmax},
PlotLegends -> Automatic, FrameLabel -> {t, x}]


Possible improvements:

• There's a typo in the definition of rhs. f[x]*NIntegrate[Abs[x - y]*f[y], {y, xmin, xmax}] should be NIntegrate[Abs[x - y]*f[x]*f[y], {y, xmin, xmax}] according to the question. – xzczd Apr 14 at 10:08
• @xzczd constant factors can be pulled out of the integral. – Roman Apr 14 at 11:03
• Oops…you're right. But, though I can't tell what the exact reason is, the numeric error in the solution seems to be a bit large? $f(x,t)=\frac{1}{10}$ is the solution, right? – xzczd Apr 14 at 11:21
• @xzczd the solution is constant in time only if you disregard the boundaries. The OP is not very clear on this topic. I've made all functions Piecewise[{{..., xmin <= x <= xmax}}] instead of allowing nonzero values outside of the domain $[x_{\text{min}},x_{\text{max}}]$. In this case, the constant function no longer solves the equation, because it is only constant on this interval and zero outside. Carification from the OP is needed if this is a problem. – Roman Apr 14 at 12:19
• Thanks a lot @Roman, Your explanation was convincing that as long as the cross-terms lie inside the bound, constant 1/10 is always the solution. Else, it'll slowly converge to the gaussian. However, I'm yet to figure out what it means physically..... – sammy Apr 14 at 15:54