The Finite Difference Method - An Introduction with Examples in Matlab
Florian Schmid, 2001
Note: The following examples illustrating the Finite Difference Method were created during a seminar in "Numerics of Partial Differential Equations". All Matlab files that have been used to generate the graphs are available for download.
Contents
1 Summary
The 5-Point-Stencil and the compact 9-Point-Stencil Finite Difference Discretization schemes are used to solve Poisson's equation. This is done using Matlab. The numerically computed convergence order is then compared to the theoretically predicted convergence order. Lexicographic numbering was used for all discretizations.
Poisson's Equation:
Boundary conditions:
Either Dirichlet boundary conditions:
or Neumann boundary conditions:
(where n denotes the outer normal vector) have been applied.
Notes on the downloadable Matlab files:
Each Matlab file uses a specific discretization method to compute a specific problem. To distinguish among the various files, a simple naming convention was established. The abbreviations used within the file names are shown in the table below.
An example for the naming convention would be: FpointDiC2.m - This function uses 5-Point-Stencil discretization. It deals with Dirichlet boundary conditions, and the exact analytic solution is given by equation C2.
Each Matlab function has two input arguments: n and plot. n denotes the number of subdivisions of the interval [0,1]. The discretization grid has then (n+1)^2 nodes. plot is an integer value that must contain either 0 or 1. If plot is set to 1, the file returns two plots. The first plot shows the exact solution, the second plot shows the numerically computed solution. If plot is set to 0, the two plots will not be shown.
The Matlab functions have an output argument as well. This is the maximal absolute difference between the computed solution and the exact solution. The maximum is taken over all the nodes of the grid.
| Abbreviation |
Meaning |
| Npoint |
Discretization with the 9-Point-Stencil |
| Fpoint |
Discretization with the 5-Point-Stencil |
| Di |
Dirichlet boundary conditions |
| Neu |
Neumann boundary conditions |
| Sym |
Symmetric Difference is used to discretize the Neumann boundary |
| 1 |
The exact solution of the computed problem is
|
| C2 |
The exact solution of the computed problem is
|
| C3 |
The exact solution of the computed problem is
|
| C4 |
The exact solution of the computed problem is
|
|
2 The Matlab .m Files
The following table lists the available Matlab files. Please click on the file name to download. The functions are commented extensively, therefore, no detailed explanations are given here.
| File Name |
Stencil Used |
Exact Solution |
Applied Boundary Conditions |
Remarks |
| FpointDi1.m |
5-Point |
Function 1 |
Dirichlet |
- |
| FpointDiC2.m |
5-Point |
Function C2 |
Dirichlet |
- |
| NpointDi1.m |
9-Point |
Function 1 |
Dirichlet |
Modified right side |
| NpointDiC2.m |
9-Point |
Function C2 |
Dirichlet |
Modified right side |
| NpointDiC3.m |
9-Point |
Function C3 |
Dirichlet |
Modified right side |
| NpointDiC4.m |
9-Point |
Function C4 |
Dirichlet |
Modified right side |
| FpointNeu1.m |
5-Point |
Function 1 |
Neumann |
Discretization of the Neumann boundary condition with the forward/backward difference quotient |
| FpointNeuSym1.m |
5-Point |
Function 1 |
Neumann |
Discretization of the Neumann boundary condition with the symmetric difference quotient |
| FpointNeuSymC2.m |
5-Point |
Function C2 |
Neumann |
Discretization of the Neumann boundary condition with the symmetric difference quotient |
|
3 Determination of the Convergence Order
Another Matlab file called order.m is used to compute the convergence order of the finite difference methods from above.
The output of the function is a graph that shows the maximum error over the grid step width h. For the chosen double-logarithmic scale, one obtains a straight line. The slope of the straight line gives the convergence order.
For a more detailed explanation of order.m, please see the commenting of the file.
The following graphs show all resulting convergence orders for the files above (see Section 2).
As expected: The convergence order is 2
|
A small surprise here: Normally, only a convergence order of 2 is assured if the solution is smooth enough. However, the convergence theory for Finite Elements, which is equivalent here, justifies the good result.
|
The approximation of the Neumann boundary condition has only a order of 1. Thus it dominates the convergence of the 5-Point-Stencil, which should be 2 (as seen above for FpointDi1).
|
Now, the approximation of the boundary condition is better (order 2). This results in a convergence order of 2 for the whole method.
|
The approximation of the boundary condition is of order 2. This results in a convergence order of 2 for the whole method. Surprising is the fact that a function in C2 is already smooth enough to converge with an order of 2.
|
As expected: The convergence order is 4.
|
The function in C2 is not smooth enough to converge with the expected convergence order of 4. This results from a higher error at the grid points that are close to the line where the second derivate of the function is not differentiable. How close a set of grid points is lying to this line depends on the grid step width. This also explains the shaky look of the curve.
|
This plot shows the error u-u_h for a 90x90 grid. One can observe that the error is much higher close to the line where the function is not smooth.
|
Surprise: The convergence oder is already 4. This can be explained with the Finite Element convergence theory, which can be applied in this case.
|
See explanation for NpointDiC3.
|
|
|
