5) Finite Differences

Last time:

• Introduction to programming in Julia

• Errors

• Plotting

Today:

1. Introduction to Finite Difference schemes
   1.1 Introduction
   1.2 Getting Started

1. Introduction to Finite Difference schemes

1.1 Introduction

This lecture introduces finite-difference grids, notation, and the construction of difference approximations to PDEs. We also show basic boundary-condition treatments and immediately implement a first solver to build numerical intuition. Some of the material of this section is largely based on the Thomas (1998), Numerical Partial Differential Equations: Finite Difference Methods book.

1.2 Getting Started

Consider the initial–boundary value problem (IBVP) for the heat equation in 1D

(6)\[\begin{align} v_t &= \nu v_{xx}, && x\in(0,1),\ t>0, \tag{1.2.1}\\ v(x,0) &= f(x), && x\in[0,1], \tag{1.2.2}\\ v(0,t) &= a(t),\quad v(1,t)=b(t), && t\ge 0. \tag{1.2.3} \end{align}\]

where \(f(0)=a(0)\), and \(f(1)=b(0)\).

Of course, this problem can be solved analytically. For purposes of illustration, we shall solve it numerically.

Definition

• Discretization: In applied mathematics, discretization is the process of translating continuous functions, models, variables, and equations into discrete counterparts (that a computer can digest).

• Of course, some essential information can be “lost in translation”. This would lead to undesirable outcomes.

• In this course we will learn that some discretization methods (or schemes) are more suitable than others; and how the choice of the “right” method is often problem-dependent.

• Most importantly, we will learn how to assess the goodness of numerical discretizations (in terms of accuracy, reliability, robustness, efficiency, and so on)

We begin by discretizing the spatial domain by placing a grid over the domain. For convenience, we will use for now a uniform grid (i.e., with all equally distant grid points), with grid spacing \(\Delta x = 1/M\) (see figure above).

To refer to points in the grid, we shall call the points \(x_k\), \(k=0,\ldots,M\), where

\[ x_k = k \Delta x, \quad k = 0, \ldots, M . \]

Likewise, we discretize the time domain similarly by placing a uniform grid on the temporal axis with grid spacing \(\Delta t\).

The resulting grid in time-space domain is illustrated in the figure below

Notationally, we define \(u_k^n\) to be a function defined at the point \((k \Delta x, n \Delta t)\), or the lattice point \((k,n)\). \(u_k^n\) represents the approximation to the solution of the problem of interest (the heat equation above).

Numerical differentiation

How to define discrete derivatives?

From calculus, the definition of the derivative of a function at a given point.

Continuous:

\[ \frac{\partial f}{\partial x} \left( x \right) = \lim_{h \to 0} \frac{f \left( x + h \right) - f \left( x \right)}{h} \]

Discrete:

\[ \frac{\partial f}{\partial x} \left( x \right) \approx \cancel{\lim_{h \to 0}} \frac{f \left( x + h \right) - f \left( x \right)}{h} \]

Taylor series

Classical accuracy analysis assumes that functions are sufficiently smooth, meaning that derivatives exist and Taylor expansions are valid within a neighborhood. In particular,

\[ f(x+h) = f(x) + f'(x) h + f''(x) \frac{h^2}{2!} + \underbrace{f'''(x) \frac{h^3}{3!} + \dotsb}_{O(h^3)} . \]

The big-\(O\) notation is meant in the limit \(h\to 0\). Specifically, a function \(g(h) \in O(h^p)\) (sometimes written \(g(h) = O(h^p)\)) when there exists a constant \(C\) such that

\[ g(h) \le C h^p \]

for all sufficiently small \(h\).

First-order spatial derivatives

To define derivatives on our lattice grid, we use Taylor’s expansion:

\[ u_{k+1} = u(x_k + \Delta x) \approx u(x_k) + \Delta x \partial_x u \bigg\rvert_k + \frac{1}{2} \Delta x^2 \partial^2_x u \bigg\rvert_k + \frac{1}{6} \Delta x^3 \partial^3_x u \bigg\rvert_k + O(\Delta x^4) \]

Similarly,

\[ u_{k-1} = u(x_k - \Delta x) \approx u(x_k) - \Delta x \partial_x u \bigg\rvert_k + \frac{1}{2} \Delta x^2 \partial^2_x u \bigg\rvert_k - \frac{1}{6} \Delta x^3 \partial^3_x u \bigg\rvert_k + O(\Delta x^4) \]

We can define the first-order derivative at the point \(x_k\) using the forward difference:

\[ \frac{\partial u_k^F}{\partial x} \approx \frac{ u_{k+1} - u_{k}}{\Delta x} + O(\Delta x) \]

this is a first-order approximation.

Similarly, we can define the first-order derivative at the point \(x_k\) using the backward difference:

\[ \frac{\partial u_{k}^B}{\partial x} \approx \frac{ u_{k} - u_{k-1}}{\Delta x} + O(\Delta x) \]

this is a first-order approximation.

If we use one point to the right of \(x_k\) and one point to the left of \(x_k\) we have a centered difference approximation for the first-order derivative at the point \(x_k\) using the centered difference:

\[ \frac{\partial u_{k}^C}{\partial x} \approx \frac{ u_{k+1} - u_{k-1}}{2 \Delta x} + O(\Delta x^2) \]

this is a second-order approximation of the first-order derivative.

We note that the centered difference approximates the first derivative with respect to \(x\) more accurately than either of the one-sided differences. In fact, the size of what is left is \(O( \Delta x^2 )\) versus \(\Delta x\).

First-order temporal derivative

In the time domain:

\[ u^{n+1} = u(x^n + \Delta t) \approx u(t^n) + \Delta t \partial_t u {\bigg\rvert}^n + \frac{1}{2} \Delta t^2 \partial^2_t u {\bigg\rvert}^n + \frac{1}{6} \Delta t^3 \partial^3_t u {\bigg\rvert}^n + O(\Delta t^4) \]

Similarly,

\[ u^{n-1} = u(x^n - \Delta t) \approx u(t^n) - \Delta t \partial_t u {\bigg\rvert}^n + \frac{1}{2} \Delta t^2 \partial^2_t u {\bigg\rvert}^n - \frac{1}{6} \Delta t^3 \partial^3_t u {\bigg\rvert}^n + O(\Delta t^4) \]

The first-order time derivative at the point \(x^n\) can be defined using the forward difference:

\[ \frac{\partial {u^n}^F}{\partial t} \approx \frac{ u^{n+1} - u^{n}}{\Delta t} + O(\Delta t) \tag{1.2.4} \]

this is a first-order approximation.

Similarly, we can define the first-order time derivative at the point \(x^n\) using the backward difference:

\[ \frac{\partial {u^{n}}^B}{\partial x} \approx \frac{ u^{n} - u^{n-1}}{\Delta t} + O(\Delta t) \]

this is a first-order approximation.

If we use one point forward in time relative to \(x^n\) and one point backward in time relative to \(x^n\), we have a centered difference approximation for the first-order time derivative at the point \(x^n\) using the centered difference:

\[ \frac{\partial {u^{n}}^C}{\partial t} \approx \frac{ u^{n+1} - u^{n-1}}{2 \Delta t} + O(\Delta t^2) \]

this is a second-order approximation of the first-order time derivative.

We note that the centered difference approximates the first derivative with respect to \(t\) more accurately than either of the one-sided differences. In fact, the size of what is left is \(O( \Delta t^2 )\) versus \(\Delta t\).

Second-order spatial derivatives

To define a second-order spatial derivative, let’s repeat the same process. This time, we need to compute the derivative of a derivative. Let’s consider half step sizes. The approximated solution \(x_{k+\frac{1}{2}}^n\) is defined at the mid grid point \( (\left(k+\frac{1}{2} \right) \Delta x)\). Similarly, \(x_{k-\frac{1}{2}}^n\) is defined at the mid grid point \( (\left(k-\frac{1}{2} \right) \Delta x)\).

\[ \frac{\partial^2 u^C_k}{\partial x^2} \approx \frac{\frac{\partial u_{k+\frac{1}{2}}^F}{\partial x} -\frac{ \partial u_{k-\frac{1}{2}}^B}{\partial x} }{\Delta x} \approx \frac{ \frac{u_{k+1} - u_k}{\Delta x} - \frac{u_k - u_{k-1}}{\Delta x} }{ \Delta x} \approx \frac{u_{k-1} -2 u_k + u_{k+1}}{\Delta x^2} + O(\Delta x^2) \tag{1.2.5} \]

Now that we have defined the discrete derivatives in each dimention, we can discretize our IVP for the heat equation putting back the explicit dependence of the approximated solution on both independent variables, \(u_k^n\). Let’s use (1.2.4) and (1.2.5):

\[ \frac{ u^{n+1}_k - u^{n}_k}{\Delta t} = \nu \frac{u_{k-1}^n -2 u_k^n + u_{k+1}^n}{\Delta x^2} \quad \textrm{{(F.T.C.S.)}} \]

This scheme is abbrievated as F.T.C.S. = Forward in Time, Centered in Space. We can recast this to isolate the approximated solution at the new time step (i.e., what we want to estimate) on the left-hand side:

\[ u^{n+1}_k = u^{n}_k + \nu \frac{\Delta t}{\Delta x^2}\left(u_{k-1}^n -2 u_k^n + u_{k+1}^n \right) \tag{1.2.6} \]

Finally, the Initial Condition (I.C.) and Boundary Conditions (B.C.) can be approximated by:

\[ u^{0}_k = f(k \Delta x), \quad k = 0, \ldots, M \quad (\textrm{I.C.}) \tag{1.2.7} \]

\[ u^{n+1}_0 = a( (n+1) \Delta t), \quad n = 0, \ldots \quad (\textrm{B.C. on the left endpoint}) \tag{1.2.8} \]

and

\[ u^{n+1}_M = b( (n+1) \Delta t), \quad n = 0, \ldots (\textrm{B.C. on the right endpoint}) \tag{1.2.9} \]

We can see that eq. (1.2.7) gives \(u_k^0\), for \(k=0, \ldots M\); eq. (1.2.6) with \(n=0\) can be used to determine \(u_k^1\), for \(k=1,\ldots M-1\) (that is, the interior grid points only); and finally, eq. (1.2.8) and (1.2.9) can be used to determine \(u_0^1\) and \(u_M^1\).

Thus, eq. (1.2.6), (1.2.8) and (1.2.9) use the information at \(n=0\) to determine \(u\) at the first time step. Once \(u^1_k\), \(k=0, \ldots M\) is known at all spatial grid points, eq. (1.2.6), (1.2.8) and (1.2.9) can be used to determine \(u\) at \(n=2\). And, of course, this process can be continued to determine \(u_k^n, k+0, \ldots M\) to any desired time step \(n\).

It should be noted that it was not poissible to determine \(u^1_0\) and \(u^1_M\) by using eq. (1.2.6), since one of the subscripts \(k+1\) and \(k-1\) would reach out of bounds (less than \(0\) or greater than \(M\)) for either of the calculations. Hence, it will almost always be necessary to have some special treatment of the boundary, such as eq. (1.2.8) and (1.2.9).

• We call this FTCS scheme an explicit scheme, because we are able to solve for the variable at the \((n+1)\)-st time step explicitly.