Finite Difference Method for Solution of PDEs

It is possible to use finite differences of a smooth function $f$ to estimate the derivative of the function. Let $D_x f = \frac{d f}{d x}$ map a function to its derivative, and let $\Delta_h f(x) = f(x+h) - f(x)$ be the finite difference.

By Taylor’s theorem,

$\Delta_h f(x) = \sum_{i=1}^\infty D_x^i f(x) \frac{h^i}{i!}$

and so the first derivative can be estimated using either forward or backward difference:

$\frac{\Delta_h f(x)}{h} = \frac{f(x+h)-f(x)}{h} = D_x f(x) + O(h)$

$\frac{\Delta_{-h} f(x)}{-h} = \frac{f(x)-f(x-h)}{h} = D_x f(x) + O(h)$

Define the operators

$e_h = \frac{1}{2}\left(\Delta_h + \Delta_{-h} \right)$

$o_h = \frac{1}{2}\left(\Delta_h - \Delta_{-h} \right)$

Note that $e_h f(x)$ is an even function of $h$ and $o_h f(x)$ is an odd function of $h$ . Relating to the Taylor series,

$e_h f(x) = \sum_{i = 1}^\infty D_x^{2i} \frac{h^{2i}}{(2i)!}$

$o_h f(x) = \sum_{i = 1}^\infty D_x^{2i-1} \frac{h^{2i-1}}{(2i-1)!}$

Hence

$\frac{1}{h^2} e_h f(x) = \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} = D_x^2 f(x) + O(h^2)$

$\frac{1}{2h} o_h f(x) = \frac{f(x+h) - f(x-h)}{2h} = D_x f(x) + O(h^2)$

provides finite difference symmetric difference estimates of the first and second order derivatives of a function.

Application of Finite Difference Approximation to Partial Differential Equations.

One of the simplest partial differential equations is the 1-dimensional heat equation, described by

$\partial_t u = c\partial_x^2 u$

where operationally $\partial_q u = \frac{\partial u}{\partial q}$ is the partial derivative of $u$ with respect to variable $q$ .

In the current context, $c$ is some positive constant and $u = u(x,t)$ is a function of position $x \ge 0$ and time $t \ge 0$ . A complete specification requires specification of boundary conditions, which for simplicity will be taken to be $u(x,0) = f(x)$ . A solution will describe the evolution of $u$ over time.

Forward Difference Approximation

Using the forward difference for estimating $\partial_t u$ and the symmetric difference for estimating $\partial_x^2 u$ , the PDE becomes

$\frac{u(t+h,x) - u(t,x)}{h} + O(h) = c\frac{u(t,x+k)-2u(t,x)+u(t,x-k)}{k^2} + O(k^2)$

Finite methods require discretization of the underlying domain of definition. For simplicity, let $t_i = ih$ and $x_j = jk$ , and let $U_ij$ be the approximation of $u(x_j,t_i)$ satisfying

$\frac{U_{i+1,j} - U_{i,j}}{h} = \frac{U_{i,j+1} - 2U_{i,j} + U_{i,j-1}}{k^2}$

This is a finite difference equation where given the value of $U$ at time index $i$ , the value of $U$ at time index $i+1$ is explicitly given as

$U_{i+1,j} = U_{i,j} + chk^{-2} \left(U_{i,j+1} - 2U_{i,j} + U_{i,j-1} \right)$

To make the problem computable, the mesh has to be finite and so the domain of definition is a finite rectangle $[0,Nh] \times [0,Mk]$ . Hence the computation will be restricted to points $i = 0,1,\ldots,N$ and $j = 0,1,\ldots,M$ . Those points that fall outside these ranges will assumed to have value $0$ .

Under certain conditions, one of which is the requirement that $chk^{-2} \le \frac{1}{2}$ , the finite computation provides an estimate of the actual solution that converges as the mesh size of the grid approaches zero.

A simple double loop program can generate the mesh solution over the rectangle with a single loop initializing the boundary condition.

Navigation

About

Raedwulf ….