2
$\begingroup$

Let us consider the following equation : $$\partial_t u = \Delta_x u + K*u$$ where $K$ is a smooth kernel, $u(x,t)$ is the unknown and $x$ is in $\Omega$ a domain in 1d or 2d.

I want to numerically solve it by the method of lines : $$\dot{u_i} = \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} + (K * u)_i$$

Is this a good way to go ? And is there any trick for the integral part $(K * u)_i$ ? I find it very time consuming.

EDIT :

The equation is actually $$\partial_t u = \Delta_x u + u (K*u)$$ The method of lines gives $$\dot{u_i} = \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} + u_i (K * u)_i$$

$\endgroup$
2
$\begingroup$

Do a Fourier transform in space and solve the ODEs that describe the first $N$ Fourier coefficient. This is possible because the Laplacian turns into a multiplication with the wave number upon Fourier transform, and the convolution turns into a straight multiplication. This makes all operations reasonably trivial, particularly if all you have is a 1d domain.

$\endgroup$
3
  • $\begingroup$ My mistake. The reaction part is more elaborate, even if the most complicated part still is $K * u$. I'll edit my original post. $\endgroup$
    – bela83
    May 6 '15 at 8:58
  • $\begingroup$ Then the question becomes: how many grid points does your kernel span? $\endgroup$ May 6 '15 at 12:34
  • $\begingroup$ All of them. The kernel is not compactly supported, although rapidly decreasing farther from the origin. $\endgroup$
    – bela83
    May 6 '15 at 13:02
2
$\begingroup$

If your discretization of the reaction part of the problem is explicit, you can still use Prof. Bangerth's suggestion. You can quickly compute $K*u$ using the fast Fourier transform on both $K$ and $u$, multiplying the transforms and then transforming back. Then you can easily compute $u\cdot K*u$ at any grid point. However, you may run into stability issues with an explicit time discretization and a reasonably large time-step; reaction-diffusion equations are often stiff. We can probably give you more suggestions if you write out what the kernel is.

Since the kernel is rapidly decreasing far from the origin, you could also try to approximate the convolution $K * u$ using something like the fast multipole method. This idea was originally used to turn the $O(N^2)$ computation of gravitational $N$-body potentials into an $O(N\log N)$ computation of an approximate potential, but it's been extended to more general problems. FMM often shows up on lists of the top 10 greatest algorithms, so it's also worth knowing about on principle.

The fast multipole method and more general hierarchical matrix methods are often used for the numerical solution of integral equations.

$\endgroup$
1
  • $\begingroup$ The kernel is radial : $K({\bf x})=K(|{\bf x}|)$. It can be either a gaussian or a difference of two gaussians with different widths. $\endgroup$
    – bela83
    May 6 '15 at 19:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.