4
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

You should never use explicit method for the diffusion equation. Implicit is unconditionally stable and just as easy to implement. Also if you use an implicit method (like backward Euler or Crank-Nicolson) it will not matter how small d is. In fact you could use a dirac delta function if you wanted.

As far as speed goes, doing a implicit method with a tridiagonal solver (O(n) complexity) this should be very fast especially for 1D. I can send an easy Matlab script if necessary.

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ Your point about stability is important; anyone using explicit methods for equations with diffusive terms would do well to keep the stability limits on time-stepping in mind. I hesitate to say "never use explicit methods" because there can be reasonable uses of explicit methods on equations with diffusive terms (for example, a complicated model that is advection-dominated, and you only have access to the right-hand side). Implicit methods are certainly preferred. $\endgroup$
    – Geoff Oxberry
    Apr 8, 2014 at 22:35
  • $\begingroup$ 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. $\endgroup$
    – Echows
    Apr 9, 2014 at 11:52
  • $\begingroup$ @Echows What I mean is that your initial condition can be zero everywhere except at a single point. For example in Matlab notation you could write u=zeros(1,N); u(N/2)=1. In this case d goes from u(N/2) to u(N/2) on your grid; in other words its "length" is zero on the grid. Now its length isnt really zero probably, but must be smaller then the grid spacing so that it appears to have zero length on the grid. $\endgroup$
    – James
    Apr 9, 2014 at 14:58
  • $\begingroup$ @Echows ...continuing. I am a little confused at what you are trying to do because the 1D heat equation can be solved very very fast for large number of grid points. If your real d is really that small then your grid spacing dx will have to be smaller of course but that is a question of resolution. How many grid points do you need to resolve your d? How small is d? At some point Matlab may not be fast enough for your problem, but I am skeptical that this is the case for the 1D heat equation. $\endgroup$
    – James
    Apr 9, 2014 at 15:17
  • $\begingroup$ @Echows...One final comment. Have you profiled your code in Matlab? Are you initializing your arrays as sparse? If your arrays are dense then this will be a problem for large number of points!!!! $\endgroup$
    – James
    Apr 9, 2014 at 15:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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