# Improving the time integration of implicit discretized PDE with a non-linear source term

This might be a naive question, but when applying a implicit discretization to a PDE with a source term, should the source be averaged in time? For example if we take the diffusion equation with a non-linear source term,

$$u_t = u_{xx}+s(x,t,u)$$

We can apply the following central difference implicit scheme to the differential term,

$$\frac{u_j^{n+1} - u_j^n}{\Delta t} = \left[ (1-\theta) (u_{j-1} - 2u_{j} + u_{j+1}) + \theta (u_{j-1} - 2u_{j} + u_{j+1})\right] + s(x,t,u)$$

but how should $s(x,t,u)$ be treated? Should we simply take the value the $n$ time point (this is what I have always done in the past),

$$s(x,t) = s_j^n$$

or averaged over time,

$$s(x,t) = (1-\theta)s_j^{n+1} + \theta s_j^{n}$$

I am not sure it is possible to implement a time average in this way because, in general, the $n+1$ points in time are unknowns!

Is this a silly question? Or is there some way of improving the time integration of the above equation by taking averages in time?

If the right hand side were independent of $u$ then one would generally use the averaged form $$(1-\theta)s_(x,t^{n+1}) + \theta s(x,t^n).$$ In the nonlinear case you can't do that easily, as you note, but you can at least use some kind of extrapolation, for example approximate $$(1-\theta)s(x,t^{n+1},u^{n+1}) + \theta s(x,t^n,u^n) \approx (1-\theta)s\left(x,t^{n+1},u^{n}+\Delta t\frac{u^n-u^{n-1}}{\Delta t}\right) + \theta s(x,t^n,u^n) .$$ You will find more tricks like this in the ODE literature -- read up on semi-implicit methods for ODEs of the kind $\dot x(t) = A x(t) + f(x(t),t)$.

• I second that. Here's a classic paper describing IMEX methods:epubs.siam.org/doi/abs/10.1137/0732037 – GradGuy Jun 1 '13 at 5:55
• Thanks, will grab that next time I'm on campus. Knowing the names for the various techniques helps. When I started this equation I didn't know I was asking an IMEX equation. – boyfarrell Jun 2 '13 at 11:03

The most robust way is solving it implicitly, otherwise for stiff nonlinearity in the source function $s(u)$ you will have to use very small time step to maintain numerical stability.

Consider

$u_t = s(t,u)$

Then

$u^{n+1} = u^n + dt \; s(t+dt, u^{n+1})$

Here one can use average between n and n+1 in the RHS to increase the order, or use more temporal points for even higher order in time. But to make it simple let's assume first order implicit time-stepping.

Now the equation can be cast in the form

${F(u^{n+1})}=0$

which can be solved by the Newton method, and for the discretized system doing linear solves with the Krylov subspace is the best known approach, so this is the Newton-Krylov method.