hi all,

Many people are solving the one-dimensional heat (diffusion) equation
using central finite differences in space and forward Euler in time.
Then, the time step dt is restricted by the well-known stability
condition dt<dx^2/2/D.

To be able to take larger time steps, some people therefore use the
implicit backward Euler scheme. Then, the time step can indeed be
chosen larger than in the explicit case.

However, as both time integration methods are first order methods, and
the error of the full numerical scheme (including spatial and temporal
discretisation errors) is

e=O(C1*dt+C2*dx^2),

I am wondering if putting much effort in using the backward Euler scheme
(which requires to solve a linear system) has much sense. Isn't the time
step dt in this case also restricted to a condition dt<O(dx^2) due to
accuracy arguments?

Or more generally stated: does it only make sense to use implicit time
integration schemes to alleviate the stability restriction if the order
of those schemes is at least 2 (still for the 1D diffusion equation of
course)?

Crank-Nicholson is no harder to program
than the implicit forward Euler.

It has the same stability and is second-order
in time rather than first-order.

With one spatial dimension the resulting matrix
is tridiagonal, so the inversion is easy.

Olin Perry Norton

Chris wrote:


hi all,

Many people are solving the one-dimensional heat (diffusion) equation
using central finite differences in space and forward Euler in time.
Then, the time step dt is restricted by the well-known stability
condition dt<dx^2/2/D.

To be able to take larger time steps, some people therefore use the
implicit backward Euler scheme. Then, the time step can indeed be
chosen larger than in the explicit case.

However, as both time integration methods are first order methods, and
the error of the full numerical scheme (including spatial and temporal
discretisation errors) is

e=O(C1*dt+C2*dx^2),

I am wondering if putting much effort in using the backward Euler
scheme (which requires to solve a linear system) has much sense. Isn't
the time step dt in this case also restricted to a condition
dt<O(dx^2) due to accuracy arguments?


Yes and No. If you have 'rough' initial values (at time "t=0") AND
you want to follow the initial, short transient phase, you're right,
you have to take small time steps. (Note, the same situation
can occur for "positive time" if you have time dependent source
terms our the equation is non-linear)

Note, even if the solution is (looks) smooth the time step of an
explicit scheme is restricted by the highest possible frequencies
even if these aren't present in the current solution. If the time
step is too large, even tiny high frequency parts (e.g. due to rounding
errors) blow up immediately.
On the contrary, with an implicit scheme which is A-stable (or L-stable
if the eigenvalues of the operator in space are real) high frequency
components will be damped out immediately for all time steps.


Or more generally stated: does it only make sense to use implicit time
integration schemes to alleviate the stability restriction if the
order of those schemes is at least 2 (still for the 1D diffusion
equation of course)?


No, see above

Helmut Jarausch

Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany

Dear Helmut,

thanks for your answer, it is quite in line with what I already
suspected. However, only one thing is not yet fully clear to me. If I am
correct, for the specific case of the diffusion equation, the global
error of both the forward and backward euler scheme is

e=O( C1 * d^2u/dt^2 * dt - C2 * d^4u/dx^4 * dx^2 )

Using the PDE itself, this can be rewritten as

e=O( (C1*dt+C2*dx^2) * d^4u/dx^4).

As far as I know, the constants C1 and C2 are now independent of the
solution itself (but they depend on, e.g., the end time of the
integration).

If the above is correct, it is in line with what you are saying:

1)if we are integrating over _short times_, C1 will be of the same
magnitude of C2, and dt must be of the same magnitude as dx^2.

2)if we want to integrate over _longer times_ (to find, e.g., a state
near a steady-state), then I can understand that the error due to the
time integration becomes smaller and smaller, until we only have the
spatial discretisation error left in the steady-state. My intuitive
explanation for this is that if the scheme is stable, then the errors
also go to zero as time goes to infinity (or is this too simplistic?).
Hence in that case the constant C1 would probably be substantially
smaller than C2, and we may take larger time steps (with the implicit
scheme).


The only thing that I do not understand is why you are speaking about
"rough" initial values: the roughness (or derivatives) appear in both
the time and space discretization (for my specific model problem), and
according to the formula above they should hence not matter (right?).
So the roughness or derivatives determine the actual size of the error
e, but not how dt whould be chosen compared to dx.


Do you agree with my reformulation of your answer? Or is there still
something wrong with my reasoning?

Peter Spellucci wrote:


this is correct ... but it is not L-stable and may give (bounded but
unphysical) oscillations for large time steps. better: BDF2.
of course you need one additional first line, since this is a two step scheme
hth
peter

Could you briefly explain, or provide a pointer to
an explanation of, L-stable and what it means?

I've tried Googling the term, and I find lots of
material on this method or that method is or is
not L-stable or A-stable, but so far I haven't found
a good definition of these terms.

Is L short for Lyapunov?

Olin Perry Norton

this is correct ... but it is not L-stable and may give (bounded but
unphysical) oscillations for large time steps. better: BDF2.
of course you need one additional first line, since this is a two step scheme
hth
peter

Peter Spellucci wrote:


this is correct ... but it is not L-stable and may give (bounded but
unphysical) oscillations for large time steps. better: BDF2.
of course you need one additional first line, since this is a two step
scheme
hth
peter

Could you briefly explain, or provide a pointer to
an explanation of, L-stable and what it means?

I've tried Googling the term, and I find lots of
material on this method or that method is or is
not L-stable or A-stable, but so far I haven't found
a good definition of these terms.

Is L short for Lyapunov?

hi all,

Many people are solving the one-dimensional heat (diffusion) equation
using central finite differences in space and forward Euler in time.
Then, the time step dt is restricted by the well-known stability
condition dt<dx^2/2/D.

To be able to take larger time steps, some people therefore use the
implicit backward Euler scheme. Then, the time step can indeed be
chosen larger than in the explicit case.

However, as both time integration methods are first order methods, and
the error of the full numerical scheme (including spatial and temporal
discretisation errors) is

e=O(C1*dt+C2*dx^2),

I am wondering if putting much effort in using the backward Euler scheme
(which requires to solve a linear system) has much sense. Isn't the time
step dt in this case also restricted to a condition dt<O(dx^2) due to
accuracy arguments?

Or more generally stated: does it only make sense to use implicit time
integration schemes to alleviate the stability restriction if the order
of those schemes is at least 2 (still for the 1D diffusion equation of
course)?


Any comments on my simplistic vision of numerical algebra are welcome!

Chris