Nonlinear integral equation on the half line

Question

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?

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$).