Convergence order of three point central finite difference scheme

lmalenica · November 8, 2017, 08:13

For example, when we solve simple 1D Poisson equation by finite difference method, why central difference scheme on uniform grid of form:

$\frac{\partial^2 u_i}{\partial x^2}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}$

is second order method for solution convergence?

I understand why approximation of first derivative is second order (and that second derivative is also second order because cancelation of first order truncation error term on uniform grid), but I don’t understand why solution also converges with second order.

Basically, this approximation is equivalent to fitting local quadratic (p=2) polynomial through three points which should theoretically yield third (p+1) order method?

FMDenaro · November 8, 2017, 09:30

Quote:

Originally Posted by lmalenica

For example, when we solve simple 1D Poisson equation by finite difference method, why central difference scheme on uniform grid of form:

$\frac{\partial^2 u_i}{\partial x^2}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}$

is second order method for solution convergence?

I understand why approximation of first derivative is second order (and that second derivative is also second order because cancelation of first order truncation error term on uniform grid), but I don’t understand why solution also converges with second order.

Basically, this approximation is equivalent to fitting local quadratic (p=2) polynomial through three points which should theoretically yield third (p+1) order method?

It is very simple to show the order of magnitude of the local truncation error, use two Taylor expansion around i (+h and -h), expressed up to fourth order terms

lmalenica · November 9, 2017, 03:56

Thx FMDenaro, but let me explain in more details what is confusing me.

I used example of linear 1D Poisson equation (smooth exact solution) and compared convergence of few numerical methods with different order of basis (interpolation) functions. I used B-spline functions as basis functions. Order of convergence was determined for each method and each basis function by plotting graph "dx - error norm" in log-log scale and slope of the curve determines convergence order of method.

Theoretically convergence order depends of polynomial order of basis function (Taylor series expansion), and for Galerkin method (classical FEM) there is full convergence order, i.e. if "p" is order of basis function than solution converges with order "p+1", while first derivative converges with order "p" and second derivative with order "p-1". For collocation method that is not case and there is reduction in convergence order. Now I'm trying to figure way.

So I used also Finite difference method (FDM) with 3 point stencil CDS which is equivalent to approximation of second derivative by local polynomial of order p=2 (i.e. quadratic) and solution converges with order 2. By looking into Taylor expansion it is clear that first derivative approximation converges with order 2, and that second derivative converge also with order 2 (for uniform/symmetric grid and even-order polynomials as in this case).

Now I'm looking for more detailed explanation why central difference scheme on three point stencil has same convergence order for solution as for first two derivatives.

FMDenaro · November 9, 2017, 04:26

First, you are right when considering that a Lagrangian polynomial over three points is of degree 2 and therefore the accuracy of the interpolation of the function is O(h^3). A second order derivative is consequently a constant.

However, it is known that the FD second derivative has O(h^2) as it exploits the symmetry of the stencil that cancels some terms. The same as happens for the first derivative when apparently you have only 2 nodes and, hence the polynomial should be linear. That is like a fourth degree polynomial is used (5-points stencil, using also i+1/2 and i-1/2) for the second derivative.

You can generalize this reasoning as shown here (sorry it is in italian)
https://www.researchgate.net/publica...0686_Capitolo3

lmalenica · November 9, 2017, 18:27

And how can I derive formula for solution (not solutions derivatives) convergence order? And it still confuses me why is of order 2 when we are apparently using polynomial of degree 2 for approximation (so I would expect order 3).

FMDenaro · November 10, 2017, 04:00

Quote:

Originally Posted by lmalenica

And how can I derive formula for solution (not solutions derivatives) convergence order? And it still confuses me why is of order 2 when we are apparently using polynomial of degree 2 for approximation (so I would expect order 3).

Sorry, your question is not clear to me. Try yourself to develop the local truncation espression for your problem.

arungovindneelan · November 10, 2017, 06:10

Quote:

Basically, this approximation is equivalent to fitting local quadratic (p=2) polynomial through three points which should theoretically yield third (p+1) order method?

You could give a reference for this.
In my knowledge and test I have carried out quadratic polynomial gives third-order convergence only for interpolation. For the first derivative approximation, it will be (p) and for the second derivative, it will be (p-1) except symmetrical stencil where we can get order (p). The proof is straightforward from Taylors series expansion, assuming there is no stiffness in the solution that is the magnitude of higher derivate vanishes with the (grid size)^n.

lmalenica · November 10, 2017, 07:44

Quote:

Originally Posted by FMDenaro

Sorry, your question is not clear to me. Try yourself to develop the local truncation espression for your problem.

Ok, I know that I can write two Taylor expansion around $x_i+\Delta x$ and $x_i-\Delta x$ and sum them to get:

$\frac{\partial^2 u_i}{\partial x^2}=\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2}+O(\Delta x^2)$

If we substitute this expression for example into Laplace equation we get:

$\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2}+O(\Delta x^2)=0$

This shows truncation error for approximation of Laplace equation by central finite difference operator. Now I'm trying to figure out how from here express convergence order of solution (i.e. discretization error). Can we simply express solution in central point as:

$u_i=\frac{1}{2}\left(u_{i+1}+u_{i-1}\right)+O(\Delta x^4)$

From here it seams that linear interpolation from two neighbor values should be 4th order accurate for this problem when we refine grid? However it is only second order accurate. What am I missing?

Quote:

Originally Posted by arungovindneelan

You could give a reference for this.
In my knowledge and test I have carried out quadratic polynomial gives third-order convergence only for interpolation. For the first derivative approximation, it will be (p) and for the second derivative, it will be (p-1) except symmetrical stencil where we can get order (p). The proof is straightforward from Taylors series expansion, assuming there is no stiffness in the solution that is the magnitude of higher derivate vanishes with the (grid size)^n.

As I wrote before, I have also did some test by using different order B-spline functions. For interpolation of function we have convergence order (p+1-d), where (p) is polynomial order and (d) is derivative order that we approximate. For example for cubic B3 spline we get convergence order (4/3/2) for (solution/1st derivative/2nd derivative).

I have also solved 1D BVP governed by Poisson equation (also with cubic B-spline) and perform convergence analysis. In that case Galerkin formulation gives optimal convergence order (4/3/2), while collocation formulation gives convergence order (2/2/2). I found same results in the literature. Now I'm trying to figure why, and I decided to look what happens with FD method, and I'm not sure how to express convergences order of solution of BVP, and not convergence order of derivative approximation.

arungovindneelan · November 10, 2017, 08:35

Quote:

Originally Posted by lmalenica

$u_i=\frac{1}{2}\left(u_{i+1}+u_{i-1}\right)+O(\Delta x^4)$

From here it seems that linear interpolation from two neighbour values should be 4th order accurate for this problem when we refine grid? However, it is only second-order accurate. What am I missing?
.

If u expand this expression in Taylors series you will get only second-order accurate. second order derivative term survive so u will be getting second order convergence.

lmalenica · November 10, 2017, 09:15

Quote:

Originally Posted by arungovindneelan

If u expand this expression in Taylors series you will get only second-order accurate. second order derivative term survive so u will be getting second order convergence.

But we got this term from Taylor series, so it would mean that we are expanding already expanded expression and we would get:

$u_i=\frac{1}{2}\left(u_{i+1}+u_{i-1}\right)+O(\Delta x^4)=u_i+O(\Delta x^2)$

which doesn't have sense?

FMDenaro · November 10, 2017, 09:29

Be careful, you are confusing the "discretization error" from the "local truncation error".

Then, if yu consider the linear sistem A.x-q=0 and substitute an approximate solution you get A.xa-q=raand hence

A.ea=-ra

as you can see only if ||A|| =O(1) you get that the error is of the same magnitude order of the residual.

FMDenaro · November 10, 2017, 09:30

Quote:

Originally Posted by lmalenica

But we got this term from Taylor series, so it would mean that we are expanding already expanded expression and we would get:

$u_i=\frac{1}{2}\left(u_{i+1}+u_{i-1}\right)+O(\Delta x^4)=u_i+O(\Delta x^2)$

which doesn't have sense?

That has no sense

lmalenica · November 10, 2017, 11:39

Quote:

Originally Posted by FMDenaro

Be careful, you are confusing the "discretization error" from the "local truncation error".

Then, if yu consider the linear sistem A.x-q=0 and substitute an approximate solution you get A.xa-q=raand hence

A.ea=-ra

as you can see only if ||A|| =O(1) you get that the error is of the same magnitude order of the residual.

I assume that here ea is discretization error, and ra is residual because of truncation error? And how can I express convergence order of discretization error if I know order of highest truncated term?

FMDenaro · November 10, 2017, 11:49

Quote:

Originally Posted by lmalenica

I assume that here ea is discretization error, and ra is residual because of truncation error? And how can I express convergence order of discretization error if I know order of highest truncated term?

The key is to understand that you are discretizing the Laplace equation using second order approximation for the differential operators that is (2D problem with constant mesh size h)

d2fe/dx^2+d2fe/dy^2-[-4*fe(i,j)+fe(i+1,j)+fe(i-1,j)+fe(i,j+1)+fe(i,j-1)]/h^2= LTE of O(h^2)

where fe=fe(x,y) is the exact solution (I simplified fe(xi,yj) as fe(i,j)).

If you have that fe(x,y) is a known analytical function in your test, you conversely check the convergence of (fn is the vector of the numerical solution:

(fe -fn)= discretization error

selig5576 · November 10, 2017, 11:50

Maybe I am not understanding the question, but let me try to answer. If you have a second order approximation for u'' then your solution will typically converge on the order of 2. In steady state, such as solving u'' = f, you only worry about your discretization method accuracy. In a transient problem you have to worry about your scheme being consistently second order in space and time. For example, BTCS is second order in space but only first order in time. Crank-Nicholson is second order in time and in space.

There are situations where, despite having second order discrete operators, you might only give a first order solution. For example, if you have discontinuous data, this can affect the accuracy of your solution.

I would look at the book of Leveque for more on error analysis of discrete differential operators.

lmalenica · November 10, 2017, 13:03

Quote:

Originally Posted by FMDenaro

If you have that fe(x,y) is a known analytical function in your test, you conversely check the convergence of (fn is the vector of the numerical solution:

(fe -fn)= discretization error

That is exactly what I did. I tested simple 1D Poisson equation for which I know analytical solution (and solution is smooth). I calculate some norm of discretization error for multiple $\Delta x$ and then calculate convergence order.

If we use classical Galerkin FEM and quadratic Lagrangian polynomials we will get that solution converges by order 3. If we differentiate numerical solution and compare with exact first derivative of solution we will get that 1st derivative converges by order 2, and 2nd derivative will converge with order 1. It means that in Galerkin FEM convergence order depends directly on order of Lagrangian polynomial (I am assuming ideal case when solution is sufficiently regular).

Now I'm trying to relate this with convergence of FDM. Apparently we are dealing with quadratic polynomial (or even of degree 4; because of symmetry for second derivative approximation) but our FDM is only 2nd order method. I can't figure why. I know that collocation method for same basis function will have reduced order when compared with Galerkin.

Is there some order reduction in FDM because its nature of approximation?

Quote:

Originally Posted by selig5576

If you have a second order approximation for u'' then your solution will typically converge on the order of 2.

Why?

Quote:

Originally Posted by selig5576

There are situations where, despite having second order discrete operators, you might only give a first order solution. For example, if you have discontinuous data, this can affect the accuracy of your solution.

I would look at the book of Leveque for more on error analysis of discrete differential operators.

I'm now consider only spatial discretization (for steady state elliptical problem with regular solution).

FMDenaro · November 10, 2017, 13:18

Be careful on the fact that in variational formulation you introduce the integration and thus one derivative disappears. Therefore, your shape function is never used for a second derivative.

selig5576 · November 10, 2017, 13:37

If you have a second order discrete operator, especially in the case of FDM, your solution accuracy depends on the accuracy of your discretization process. To convince yourself, do a log-log plot in the L2 error (for example) of a variety of resolutions.

arungovindneelan · November 10, 2017, 20:06

Quote:

Originally Posted by lmalenica

But we got this term from Taylor series, so it would mean that we are expanding already expanded expression and we would get:

$u_i=\frac{1}{2}\left(u_{i+1}+u_{i-1}\right)+O(\Delta x^4)=u_i+O(\Delta x^2)$

which doesn't have sense?

Now I get the problem. Order of accuracy is defined by order of the PDE (or ODE) we discretize not by the outcome of PDE ( [math]u_i[math] ).

$u_i=\frac{1}{2}\left(u_{i+1}+u_{i-1}\right)=u_i+O(\Delta x^2)$

In this expression, you could think that while expanding in Taylors series second order term survive, but Laplace equation says the second derivative is zero so it could give fourth order but unluckily numerical method approximate the solution. Numerically we didn't ensure $\frac{d^2u}{dx^2} =0$ so it will not give fourth order accuracy. A few years ago I tried this in one assignment at that time I am also got second-order accuracy.

Anyway, I like this way of thinking about this simple problem!

lmalenica · November 10, 2017, 21:26

My opinion is that polynomial order determines theoretical (maximum) order of approximation. Differential operator (i.e. order of PDE/ODE) should affect only accuracy of the approximation and not the order. From results what I see its seems that that formulation/method can reduce order of method but I don't know how to verify/prove it.

One more thing that confuses me here. If we consider unsteady convection equation and discretized both spatial and temporal derivatives with forward (or backward) finite difference our method is first order accurate both in space and time. Both terms are derived by using first two terms from Taylor series (linear approximation) which implies second order method when we consider solution (not derivative) convergence order. Moreover we know that Euler temporal discretization is locally (if we consider single time step) 2nd order accurate, and that loses one order because accumulation of error during time (multiple time steps). Why then is spatial discretization only 1st order and not 2nd order accurate?

November 8, 2017, 08:13	Convergence order of three point central finite difference scheme	#1
lmalenica New Member Luka Malenica Join Date: Feb 2017 Location: Split, Croatia Posts: 8 Rep Power: 9	For example, when we solve simple 1D Poisson equation by finite difference method, why central difference scheme on uniform grid of form: $\frac{\partial^2 u_i}{\partial x^2}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}$ is second order method for solution convergence? I understand why approximation of first derivative is second order (and that second derivative is also second order because cancelation of first order truncation error term on uniform grid), but I don’t understand why solution also converges with second order. Basically, this approximation is equivalent to fitting local quadratic (p=2) polynomial through three points which should theoretically yield third (p+1) order method? arungovindneelan likes this.

November 9, 2017, 04:26		#4
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	First, you are right when considering that a Lagrangian polynomial over three points is of degree 2 and therefore the accuracy of the interpolation of the function is O(h^3). A second order derivative is consequently a constant. However, it is known that the FD second derivative has O(h^2) as it exploits the symmetry of the stencil that cancels some terms. The same as happens for the first derivative when apparently you have only 2 nodes and, hence the polynomial should be linear. That is like a fourth degree polynomial is used (5-points stencil, using also i+1/2 and i-1/2) for the second derivative. You can generalize this reasoning as shown here (sorry it is in italian) https://www.researchgate.net/publica...0686_Capitolo3 selig5576 and lmalenica like this.

November 10, 2017, 11:50	Second order FDM	#15
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 10	Maybe I am not understanding the question, but let me try to answer. If you have a second order approximation for u'' then your solution will typically converge on the order of 2. In steady state, such as solving u'' = f, you only worry about your discretization method accuracy. In a transient problem you have to worry about your scheme being consistently second order in space and time. For example, BTCS is second order in space but only first order in time. Crank-Nicholson is second order in time and in space. There are situations where, despite having second order discrete operators, you might only give a first order solution. For example, if you have discontinuous data, this can affect the accuracy of your solution. I would look at the book of Leveque for more on error analysis of discrete differential operators.

November 10, 2017, 13:18		#17
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	Be careful on the fact that in variational formulation you introduce the integration and thus one derivative disappears. Therefore, your shape function is never used for a second derivative. lmalenica likes this.

November 10, 2017, 13:37	Convergence	#18
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 10	If you have a second order discrete operator, especially in the case of FDM, your solution accuracy depends on the accuracy of your discretization process. To convince yourself, do a log-log plot in the L2 error (for example) of a variety of resolutions.

November 9, 2017, 03:56		#3
lmalenica New Member Luka Malenica Join Date: Feb 2017 Location: Split, Croatia Posts: 8 Rep Power: 9	Thx FMDenaro, but let me explain in more details what is confusing me. I used example of linear 1D Poisson equation (smooth exact solution) and compared convergence of few numerical methods with different order of basis (interpolation) functions. I used B-spline functions as basis functions. Order of convergence was determined for each method and each basis function by plotting graph "dx - error norm" in log-log scale and slope of the curve determines convergence order of method. Theoretically convergence order depends of polynomial order of basis function (Taylor series expansion), and for Galerkin method (classical FEM) there is full convergence order, i.e. if "p" is order of basis function than solution converges with order "p+1", while first derivative converges with order "p" and second derivative with order "p-1". For collocation method that is not case and there is reduction in convergence order. Now I'm trying to figure way. So I used also Finite difference method (FDM) with 3 point stencil CDS which is equivalent to approximation of second derivative by local polynomial of order p=2 (i.e. quadratic) and solution converges with order 2. By looking into Taylor expansion it is clear that first derivative approximation converges with order 2, and that second derivative converge also with order 2 (for uniform/symmetric grid and even-order polynomials as in this case). Now I'm looking for more detailed explanation why central difference scheme on three point stencil has same convergence order for solution as for first two derivatives.

November 9, 2017, 18:27		#5
lmalenica New Member Luka Malenica Join Date: Feb 2017 Location: Split, Croatia Posts: 8 Rep Power: 9	And how can I derive formula for solution (not solutions derivatives) convergence order? And it still confuses me why is of order 2 when we are apparently using polynomial of degree 2 for approximation (so I would expect order 3).

November 10, 2017, 09:29		#11
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	Be careful, you are confusing the "discretization error" from the "local truncation error". Then, if yu consider the linear sistem A.x-q=0 and substitute an approximate solution you get A.xa-q=raand hence A.ea=-ra as you can see only if \|\|A\|\| =O(1) you get that the error is of the same magnitude order of the residual.

November 10, 2017, 21:26		#20
lmalenica New Member Luka Malenica Join Date: Feb 2017 Location: Split, Croatia Posts: 8 Rep Power: 9	My opinion is that polynomial order determines theoretical (maximum) order of approximation. Differential operator (i.e. order of PDE/ODE) should affect only accuracy of the approximation and not the order. From results what I see its seems that that formulation/method can reduce order of method but I don't know how to verify/prove it. One more thing that confuses me here. If we consider unsteady convection equation and discretized both spatial and temporal derivatives with forward (or backward) finite difference our method is first order accurate both in space and time. Both terms are derived by using first two terms from Taylor series (linear approximation) which implies second order method when we consider solution (not derivative) convergence order. Moreover we know that Euler temporal discretization is locally (if we consider single time step) 2nd order accurate, and that loses one order because accumulation of error during time (multiple time steps). Why then is spatial discretization only 1st order and not 2nd order accurate?

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Code verification of second order scheme	kmandar	SU2	0	June 26, 2016 10:28
AUSM scheme ? Central Scheme	boling	Main CFD Forum	7	January 7, 2016 02:41
High-order finite difference scheme for incompressible Navier Stokes equation	novan tofany	Main CFD Forum	2	April 5, 2012 06:34
Force can not converge	colopolo	CFX	13	October 4, 2011 22:03
Definition of limiter function for central dirrerencing scheme	sebastian_vogl	OpenFOAM Running, Solving & CFD	0	January 5, 2009 11:08