# The most efficient way to solve diffusion equation with concentrated initial condition

I want to solve the diffusion equation, i.e. $$\dot{f} - f'' = 0$$ with a boundary condition $f(0) = f(1) = 0$ and with an initial condition that $f$ is a boxcar function concentrated over some small region of size $d \ll 1$. The units are scaled so that the diffusion constant is equal to unity. Then I want to use the resulting solution to compute an integral of the form $$\int dt \int dx F(f(x), f'(x)) ,$$ where $F$ is some function (it is actually quadratic in $f$ with some $x$-dependent coefficients if this makes any difference).

I have tried two different ways to do this. First is to solve the diffusion equation analytically in Fourier space and then write my integral also in Fourier space and calculate it numerically. The second is to solve the diffusion equation numerically (I've tried both explicit and implicit methods and the well-known Crank-Nicolson method).

In all approaches, I end up with the same problem: The initial condition for $f$ is very concentrated which means I need a ridiculously small space and time discretization (corresponds to taking into account very large Fourier modes in F-space solution). On the other hand, since I need to solve the equation also for large times (to calculate the integral), I need to take into account also solutions in large times. Fourier-space method is not faster because in that I actually need to compute three integrals numerically (two Fourier sums with a kernel that is calculated from the $x$-integral). This makes my calculations very slow.

Are there any better methods to do what I just described? Since diffusion equation is very simple, I feel pretty stupid for not being able to do this faster. I've been thinking about implementing some kind of adaptive grid that would be concentrated to small times and close to the initial condition, but this seems like a complicated thing to do tho solve such a simple problem.

• What do you mean it doesn't matter how small $d$ is? I must still have $dx \ll d$ which slows my computations if I set $d$ very small. Also, for implicit method, to get a good accuracy I must have $dt \ll d^2$ (otherwise the fast relaxation close to zero-time is calculated inaccurately). These two conditions combined tell me that I must have $dt \sim dx^2$ which is, up to a constant, the stability condition for the explicit method. Since the implicit method involves the additional step of matrix inversion, it is not clear to me that it is always faster than the explicit method. Apr 9, 2014 at 11:52