Approximations used in Crank-Nicolson method for solving PDEs numerically

[pafcu]

Hi all,

A diffusion equation is of the form
\frac{\partial u}{\partial t}=k \frac{\partial^2u}{\partial x^2}

Usually an equation like this seems to be solved numerically using the Crank-Nicolson method:
\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = \frac{k}{2 (\Delta x)^2}\left((u_{i + 1}^{n + 1} - 2 u_{i}^{n + 1} + u_{i - 1}^{n + 1}) + (u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n})\right)

i.e the time-related part is approximated using a 2-point forward difference approximation and the space-part using an average of two 3-point central difference approximation. (u_{i}^{n} means value of u at point i at timestep n, i.e u(i,n))

Is there a specific reason why a central difference is not used for the time? Since a central difference results in a smaller error, yet it's not usually used, what drawback does it have?

What about the central difference for the space-related part? Is there a reason for using the 3-point version instead of the 5-point version? Or is the simpler one usually used because it results in a tridiagonal matrix which is easier to solve?

edit: Clarified notation

 

[John Creighto]

Is there a specific reason why a central difference is not used for the time? Since a central difference results in a smaller error, yet it's not usually used, what drawback does it have?

I don't know anything about this method but How would know what the value is in the future time step before the derivative was estimated? (I think though Rung Kutta Methods for ODE's try and do what you suggest).

 

[pafcu]

I don't know anything about this method but How would know what the value is in the future time step before the derivative was estimated? (I think though Rung Kutta Methods for ODE's try and do what you suggest).

I was unclear. The point is that I wish to know the value of u at the next time step at every point, solving the mentioned equation for u_{i}^{n + 1} for each i (requires solving a set of linear equations) tells me just that.

 

[John Creighto]

I was unclear. The point is that I wish to know the value of u at the next time step at every point, solving the mentioned equation for u_{i}^{n + 1} for each i (requires solving a set of linear equations) tells me just that.

In that case I'm surprised you have enough equations given it looks like there will be at least two unknowns per equation. However, I suppose you use an integration step to get that other unknown.

 

[pafcu]

In that case I'm surprised you have enough equations given it looks like there will be at least two unknowns per equation. However, I suppose you use an integration step to get that other unknown.

No, number of unknowns is not a problem.

Assuming u is defined at m m points (it has been discretized) we obviously end up with m unknowns: u_i^{n+i}, 1 \leq i \leq m.
We get m-2 equations by plugging in i=2,3,4,...,m-1 into the equation in the original post. Plugging in i=1 doesn't work because u_{0}^{n+1} required by the approximation isn't defined. Similarly you can't plug in the value i=m because u_{m+1}^{n+1} isn't defined.

Is all hope lost? No, because if we additionally specify two boundary conditions (say u_0^{n+1}=0 and u_m^{n+1}=0) we end up with the required m equations and the system has a unique solution.

There is no question this method doesn't work. My specific problem was with some of the details, why certain approximations are used instead of others that should result in more accurate results.

 

[John Creighto]

No, number of unknowns is not a problem.

Assuming u is defined at m m points (it has been discretized) we obviously end up with m unknowns: u_i^{n+i}, 1 \leq i \leq m.
We get m-2 equations by plugging in i=2,3,4,...,m-1 into the equation in the original post. Plugging in i=1 doesn't work because u_{0}^{n+1} required by the approximation isn't defined. Similarly you can't plug in the value i=m because u_{m+1}^{n+1} isn't defined.

Is all hope lost? No, because if we additionally specify two boundary conditions (say u_0^{n+1}=0 and u_m^{n+1}=0) we end up with the required m equations and the system has a unique solution.

There is no question this method doesn't work. My specific problem was with some of the details, why certain approximations are used instead of others that should result in more accurate results.

I should have realized that. Anyway, if you wanted to do the central difference for the derivative, you would also need u_{i}^{n+2} if you didn't want to introduce any time legs. This means you would need twice as many equations. This would obviously slow down your algorithm but it might help accuracy.

 

[Nick Stokes]

Yes, there is a reason. It would be unstable. A quick way to check. The instability mode is usually an alternating wave, +-+-+-+... It will be modified by a factor h say in each step. Substitute that in the algorithm - ie in successive timesteps, (+1,-1,+1...), (+h,-h,+h..), (+h^2,-h^2,+h^2..).
For a two-step algorithm like C-N, you'll get a linear equation for h, which will be <1. The wave decays, which is what you want and expect.
For a three-step algorithm like the one you propose, you'll get a quadratic in h. You'll quickly see that the product of the roots is 1. One is less than 1, and is also the one you expect. In fact, it's more accurate than C-N. But there is a spurious second root that is greater than 1. That causes such waves to grow exponentially, and that is why the solutions are unstable. It's a common problem with multi-step explicit algorithms.

 

[tanujkush]

Hi all,

A diffusion equation is of the form
\frac{\partial u}{\partial t}=k \frac{\partial^2u}{\partial x^2}

Usually an equation like this seems to be solved numerically using the Crank-Nicolson method:
\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = \frac{k}{2 (\Delta x)^2}\left((u_{i + 1}^{n + 1} - 2 u_{i}^{n + 1} + u_{i - 1}^{n + 1}) + (u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n})\right)

i.e the time-related part is approximated using a 2-point forward difference approximation and the space-part using an average of two 3-point central difference approximation. (u_{i}^{n} means value of u at point i at timestep n, i.e u(i,n))

Is there a specific reason why a central difference is not used for the time? Since a central difference results in a smaller error, yet it's not usually used, what drawback does it have?

What about the central difference for the space-related part? Is there a reason for using the 3-point version instead of the 5-point version? Or is the simpler one usually used because it results in a tridiagonal matrix which is easier to solve?

edit: Clarified notation

The reason you go backward in time yet a central difference in space is precisely due to the ease involved in solving. If you went central in time as well, you would have an extra u evaluated at n-1, which would disturb the tri-diagonal structure of the RHS.

The reason the 3 point scheme is used for space is also for the same reason that you mentioned, for ease in computation. There are problems where a 5-point difference scheme is used, but then the physics is different there.

As a good exercise, you could try using a 3 point difference scheme for time and compare your results.

Another reason I came up with is that in a diffusion problem, the state at any time t would only depend on time t-1. Think about it, the temperature at a point now does not care what the temperature at that same point was say ten delta time back. All in all, its just physics!