Finite Difference Method I: Discretization, Elliptic Problems, and Numerical Linear Algebra

Last Update: 2026-09-10

Purpose

This part is not intended to make you an expert in the finite difference method (FDM).

Its purpose is to establish the basic numerical concepts and implementation skills needed for later study of:

• time-dependent partial differential equations;
• hyperbolic and propagation-dominated problems;
• high-order methods;
• discontinuous Galerkin (DG) and Runge–Kutta DG methods;
• machine-learning-assisted numerical methods.

Finite differences are used here because their discrete structure is transparent and easy to experiment with.

The main learning objective is not to memorize formulas, but to understand the relationship

Learning Levels

Core

The following topics must be implemented and understood:

• Taylor expansion and finite-difference stencils;
• local truncation error and order of accuracy;
• consistency, stability, and convergence;
• the 1D Poisson equation;
• basic treatment of Dirichlet and Neumann boundary conditions;
• error norms and convergence-rate studies;
• basic sparse linear systems.

Optional

The following topics are recommended but are not required for progressing to DG:

• the 2D Poisson equation;
• iterative linear solvers;
• preconditioning;
• higher-order boundary formulas;
• sparse matrix storage formats;
• basic complexity and conditioning experiments.

Awareness

Conceptual understanding is sufficient:

• multigrid;
• advanced preconditioners;
• finite-volume methods;
• finite-element methods;
• DG methods;
• large-scale sparse direct solvers.

1. Finite-Difference Approximation

1.1 Uniform grid

Let

where $h$ is the grid spacing.

For a sufficiently smooth function $u(x)$, denote $u_i \approx u(x_i).$ The basic goal is to approximate derivatives using values of $u$ at neighboring grid points.

1.2 First derivatives

Taylor expansion gives

Similarly,

From these expansions:

Forward difference

Backward difference

Central difference

The central difference is formally second-order accurate, but it is not automatically suitable for every time-dependent PDE. Stability and qualitative behavior must also be considered.

1.3 Second derivatives

The standard central approximation is

Taylor expansion gives

The stencil is therefore second-order accurate.

1.4 Higher-order stencils

Higher-order finite-difference formulas can be obtained using:

• Taylor expansion;
• the method of undetermined coefficients;
• polynomial interpolation;
• Fourier-symbol analysis.

For example, a fourth-order approximation to the first derivative is \(u'(x_i) \approx \frac{ -u_{i+2}+8u_{i+1}-8u_{i-1}+u_{i-2} }{12h}.\)

You should not only derive the formula, but also verify numerically that the error decreases as $\mathcal{O}(h^{4})$ for a smooth test function.

2. Truncation Error, Consistency, Stability, and Convergence

These concepts must be kept distinct.

2.1 Local truncation error

The local truncation error is obtained by substituting the exact solution into the discrete equation.

For example, for the second derivative,

For a smooth function,

This measures the local error of the discrete operator.

2.2 Discrete or global error

The numerical solution $u_i^{h}$ is compared with the exact solution at the grid points:

Typical discrete norms include

and the grid-weighted discrete $L^{2}$ norm

The observed convergence rate between two grids is

2.3 Consistency

A discretization is consistent if its discrete operator approaches the continuous differential operator as $h\to 0$.

Consistency is a local property of the approximation.

2.4 Stability

Stability concerns the behavior of perturbations caused by:

• perturbations in the data;
• round-off errors;
• initial errors;
• linear-system solution errors;
• time-stepping errors.

A consistent method can still be useless if it is unstable.

2.5 Convergence

A method is convergent if the numerical solution approaches the exact solution as the grid is refined.

In many standard settings, consistency together with a suitable stability property leads to convergence. You should understand this relationship conceptually, without treating it as a universal theorem that applies automatically to every discretization.

3. One-Dimensional Poisson Equation

Consider \(\begin{cases} -u''(x)=f(x), & x\in(0,1),\\ u'(0)=\alpha,\\ u(1)=\beta. \end{cases}\)

This problem is used to practice:

• finite-difference discretization;
• boundary-condition treatment;
• matrix assembly;
• linear-system solution;
• convergence verification.

3.1 Manufactured solution

For verification, choose

Then

and the boundary data are

Thus the test problem is

3.2 Grid and interior discretization

Let

At an interior grid point $x_j$,

Therefore,

Equivalently,

3.3 Neumann boundary condition

A simple first-order approximation is

This is easy to implement, but the boundary approximation is only first-order accurate and may reduce the observed global convergence rate.

For a second-order treatment, use

The second-order formula is recommended for the convergence experiment.

The student should compare the two treatments and observe how boundary discretization affects the global error.

3.4 Dirichlet boundary condition

At $x=1$,

The known boundary value is incorporated into the right-hand side of the neighboring interior equation.

3.5 Linear system

After applying the boundary conditions, the discrete problem has the form

Depending on the scaling convention, $A$ can be assembled using either:

or the scaled matrix

The right-hand side must be scaled consistently.

You should understand:

• which unknowns are included in the system;
• how known boundary values modify the right-hand side;
• how the Neumann equation changes the first row;
• why the resulting matrix is sparse and banded.

4. Implementation and Verification

Core implementation

Implement the 1D Poisson solver in Python or C++.

The implementation should include:

1. grid generation;
2. evaluation of the exact solution and forcing term;
3. matrix or stencil assembly;
4. Dirichlet and Neumann boundary treatment;
5. linear-system solution;
6. error computation;
7. convergence-rate estimation;
8. visualization.

Required experiments

Ex 1: Interior discretization

Use a sufficiently accurate boundary treatment and verify the expected second-order convergence.

Ex 2: Neumann boundary treatment

Compare:

with

Explain the observed difference in global convergence.

Ex 3: Error norms

Compare:

• $L^{\infty}$ error;
• discrete $L^{2}$ error;
• optionally, a discrete $H^{1}$-type seminorm.

Ex 4: Grid refinement

Use at least four grid sizes and report:

5. Sparse Linear Algebra

This section is partly optional but provides useful computational background.

5.1 Direct solvers

For a tridiagonal system, the Thomas algorithm can solve the system in approximately $\mathcal{O}(N)$ operations.

General dense Gaussian elimination is much more expensive and should be used only for comparison on small systems.

5.2 Iterative solvers

The student should be aware of:

• Jacobi iteration;
• Gauss–Seidel iteration;
• Successive Over-Relaxation;
• Conjugate Gradient for symmetric positive-definite systems;
• GMRES for general nonsymmetric systems.

For the standard 1D Poisson problem with compatible boundary treatment, the resulting system is typically sparse and well suited to tridiagonal direct solvers or iterative methods.

The student should compare:

• residual reduction;
• error reduction;
• iteration count;
• computational time;
• dependence on grid resolution.

5.3 Conditioning

As the grid is refined, the condition number of the discrete Poisson matrix increases. The student should understand qualitatively that:

• finer grids improve spatial resolution;
• but the associated algebraic system becomes more difficult to solve;
• iterative methods may require preconditioning;
• discretization cost and linear-algebra cost are connected.

A detailed condition-number analysis is Optional.

6. Optional: Two-Dimensional Poisson Equation

Consider

For example, impose Dirichlet conditions on the horizontal boundaries and Neumann conditions on the vertical boundaries:

Be careful that the outward normal has opposite directions on the two vertical boundaries:

On a uniform grid,

the standard five-point approximation at an interior point is

For $h_x=h_y=h$,

After ordering the unknowns, the matrix has a sparse block-tridiagonal structure. For compatible treatment of the mixed boundary conditions and after eliminating prescribed Dirichlet values, the resulting system is typically sparse and, in the standard formulation, symmetric positive definite.

Optional goals

• assemble the five-point stencil;
• handle mixed boundary conditions;
• use sparse matrix storage;
• solve with a sparse direct solver or CG;
• study computational cost as $N_x,N_y$ increase.

Multigrid should be treated as an awareness topic at this stage rather than a required implementation.

7. Numerical Experiment Standards

Every computational experiment should state:

• the PDE and boundary/initial data;
• the exact or reference solution;
• the spatial grid;
• the time step, if applicable;
• the stability-related parameter;
• the error norm;
• the computational cost;
• the observed numerical behavior;
• the interpretation of the result.

You should distinguish between:

• discretization error;
• time-integration error;
• boundary-condition error;
• linear-system solution error;
• round-off error.

A plot without an explanation is not considered a complete numerical experiment.

8. Connection to Later DG and Machine Learning Work

The concepts introduced here will reappear in later work:

The purpose of this transition is not to replace FDM by DG immediately. Rather:

FDM makes the basic numerical phenomena transparent; DG will later provide a more flexible high-order framework for studying local structure, shocks, unresolved scales, and numerical post-processing.