I want to numerically solve the following equation for $\phi$ on $\mathbb{R}_+^{*2}$:

$$ \partial_t \phi (w, t) = \int_0^{+ \infty} k(\alpha w + \beta w', w') \phi(\alpha w + \beta w', t) \phi(w', t) ~dw' - \int_0^{+ \infty} k(\gamma w + \delta w', w') \phi(\gamma w + \delta w', t) \phi(w', t) ~dw' $$ $$ \phi(w, 0) = \phi_0(w) $$ $$ \partial_t \phi(0, t) = 0 $$

with $k$ a given kernel and $\alpha, \beta, \gamma, \delta$ four given real numbers. The initial condition $\phi_0$ and the kernel $k$ are assumed to have sufficient regularity.

I expect interesting solutions to decrease as $w$ grows so I will define my domain as $(0, W) \times (0, T)$.

I would like to try a semi-discretization in $w$ only, so that I can feed a system of time-evolution equations to an ODE solver. However I do not know how to treat the integrals and it makes me wonder if it is at all possible.

Is it possible, and if yes, what method could I try?

  • 1
    $\begingroup$ Should we assume that the integral operators are local in time (no time delays, etc.)? Also, is the kernel local, $k(w,w') = k(w-w')$, or must we allow it to be nonlocal in general? And do you know anything special about $\alpha,\beta,\gamma,\delta$ or must we take them as arbitrary real (and positive?) parameters? $\endgroup$
    – Endulum
    Commented Apr 20, 2022 at 12:31
  • $\begingroup$ Yes, the integral operators are local in time, I will reword the r.h.s. On the other hand, $k$ is not local. In my problem, $\alpha=\gamma=1$, $\beta<0$ and $\delta>0$. $\endgroup$
    – bela83
    Commented Apr 20, 2022 at 15:03
  • $\begingroup$ For semidiscretization, you will need to replace the integrals by quadrature. A common choice would be to use a subdivision in $w$-space into small elements, and then the use of a trapezoidal rule to approximate the integral with it. $\endgroup$ Commented Apr 20, 2022 at 23:31
  • $\begingroup$ Would you comment on whether $k$ is a singular kernel? $\endgroup$ Commented Apr 22, 2022 at 22:21
  • $\begingroup$ @Aruralreader I am unsure what it means in that context but as stated in the question, $k$ can be assumed to be - say - $\mathcal{C}^{\infty}$ on its domain. $\endgroup$
    – bela83
    Commented Apr 25, 2022 at 8:30

1 Answer 1


As mentioned in comments, the first thing to try with this problem is to discretize the integrals with a simple quadrature rule and shove the semi-discrete form into a canned ODE integrator. I don't claim that the following is an efficient method or even a good one, but it is the simplest thing that might work and should serve as a starting point for researching fancier approaches if the need arises.

The basic strategy will be to discretize the integrals with the trapezoid rule and then use linear interpolation to handle evaluating $\phi$ between grid points. To start, define $$ I_{\alpha,\beta}(\phi, w) = \int_0^\infty k(\alpha w+\beta w',w')\phi(\alpha w+\beta w')\phi(w')dw' $$ (with the time dependence suppressed). We truncate the domain from $[0,\infty)$ to $[0,W]$ and subdivide this integration region into $N-1$ equal-size panels of width $h=W/(N-1)$, applying the trapezoid rule to each $$ I_{\alpha,\beta}(\phi,w) \approx \frac{h}{2} \left[ k(\alpha w,0)\phi(\alpha w)\phi(0) + k(\alpha w + \beta W) \phi(\alpha w + \beta W) \phi(\beta W)\right] + h\sum_{j=2}^{N-1} k(\alpha w+\beta w_j, w_j)\phi(\alpha w+\beta w_j)\phi(w_j) $$ where $w_i = (i-1)h$, so that $w_1=0$ and $w_N=W$. Normally, we would like the vector $\phi_i=\phi(w_i)$ to be our unknowns, but this is tricky because our quadrature rule asks for $\phi$ to be evaluated between grid points. This happens because the approximation for $I_{\alpha,\beta}(\phi,w_i)$ depends on $\phi(\alpha w_i+\beta w_j)$. Linear interpolation is one way to address this, and it is consistent with the use of the trapezoid rule, since both are based on a piecewise linear approximation to $\phi(w)$. Most numerical software packages have ready-made linear interpolation routines.

From a semi-discrete point-of-view, your implementation of the RHS function could be as simple as evaluating the trapezoid rule and interpolating where needed. One aspect that has been glossed over so far is the boundary conditions at $w=0$ and $w=W$. You will need to make sure those are properly accounted for in evaluating the trapezoid rule. You will also need to handle the case where the interpolated point $w_*$ lies off your grid ($w_*<0$ or $w_*>W$).

  • $\begingroup$ Thanks Endulum it looks like a really good starting point. I was not sure about linear interpolation, so thanks for clarifying and justifying. BC at $w=0$ is clear for me, it will be Neumann homogeneous. BC at $w=W$ is less clear, I will have to try and see. $\endgroup$
    – bela83
    Commented Apr 21, 2022 at 8:46
  • $\begingroup$ I hope it's useful to get you started. For the BC at $W$, if your solution is effectively zero there, you might be able to get away with a homogeneous Dirichlet condition. Maybe try that and then follow up here if it doesn't work out. $\endgroup$
    – Endulum
    Commented Apr 21, 2022 at 14:38

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.