A mimimum of theory for FD (finite difference) methods

1  The problem

We consider here a special type of a two point boundary value problem: Find a function $y$ of the real variable $x$ which satisfies the ordinary differential equation

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{y"}& =\hfill & f(x,y,y\text{'})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{with}f:\hspace{0.5em}[a,b]\times R\times R\to R\hspace{0.5em},\hfill \\ \multicolumn{1}{c}{}\end{array}}$

and the boundary conditions

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{g_1(y(a),y(b),y\text{'}(a),y\text{'}(b))}& =\hfill & 0\hfill \\ \multicolumn{1}{c}{g_2(y(a),y(b),y\text{'}(a),y\text{'}(b))}& =\hfill & 0\hspace{0.5em}.\hfill \end{array}}$

Here $f,\hspace{0.5em}{g}_1,g_2$ are given functions of 3 respectively 4 variables. The simplest model for the boundary conditions is clearly

${\displaystyle y(a)\hspace{0.5em}=\hspace{0.5em}y(b)\hspace{0.5em}=\hspace{0.5em}0}$

which is also used in some of our examples. For general cases of $f,\hspace{0.5em}{g}_1,\hspace{0.5em}{g}_2$ there exist no existence results, not to speak about uniqueness. Indeed in our examples there are two quite harmless looking cases which have two solutions and cases of nonsolvability can easily be constructed. Hence the situation here is quite different from the case of initial value problems. We assume in the following that a solution exists which is also sufficiently often differentiable in the interior of $]a,b[$, typically we will need four derivatives of $y$ at least.

2  Solution by discretization

We now replace the continuous problem by one in finitely many unknowns by a technique known as discretization by finite differences:

Step 1:

We choose a grid of constant width $x_i=a+i·h$ for $i=0,\dots ,N$ and $h=\frac{b-a}{N}$. We write now

${\displaystyle y"(x_i)=f(x_i,y_i,y{\text{'}}_i)}$

with the simplifying notation $y_i:=y(x_i)$, $y{\text{'}}_i:=y\text{'}(x_i)$.

Step 2:

We replace the derivatives of $y$ at $x_i$ by suitable derivatives of interpolation polynomials of some neighboring points $(x_{i-j},y_{i-j}),...,(x_{i+k},y_{i+k})$. For reasons concerning the properties of the resulting system of equations the degree chosen usually is 1 or 2. For example with $j=k=1$

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{y_i\text{'}}& =\hfill & \frac{y_{i+1}-y_{i-1}}{2h}+O(h^2)\hfill \\ \multicolumn{1}{c}{y_i"}& =\hfill & \frac{y_{i+1}-2y_i+y_{i-1}}{h^2}+O(h^2)\hfill \end{array}}$

Step 3:

Let $y_i^h$ denote the approximation for $y(x_i)$. We get, by neglecting the $O-$terms the following system defining it:

${\displaystyle \oplus \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\frac{y_{i+1}^h-2y_i^h+y_{i-1}^h}{h^2}-f(x_i,y_i^h,\frac{y_{i+1}^h-y_{i-1}^h}{2h})\hspace{0.5em}=\hspace{0.5em}0\hspace{0.5em},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}1\le i\le N-1}$

Step 4:

The treatment of the boundary conditions. We first stick on the homogeneous conditions

${\displaystyle y(a)\hspace{0.5em}=\hspace{0.5em}y(b)\hspace{0.5em}=\hspace{0.5em}0}$

resulting in

${\displaystyle y_0^h\hspace{0.5em}=\hspace{0.5em}{y}_N^h\hspace{0.5em}=\hspace{0.5em}0\hspace{0.5em}.}$

This leads to a system of $N+1$ coupled equations for the unknown approximations $y_i^h$ for the values $y(x_i)$. This system will be nonlinear if the function $f$ is nonlinear in $y$ or $y\text{'}$. Example $N=4$, i.e. $i=0,1,2,3,4$:

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{i=1:\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\frac{y_2^h-2y_1^h+y_0^h}{h^2}}& =\hfill & f(x_1,y_1^h,\frac{y_2^h-y_0^h}{2h})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{with}{y}_0^h=0\hfill \\ \multicolumn{1}{c}{i=2:\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\frac{y_3^h-2y_2^h+y_1^h}{h^2}}& =\hfill & f(x_2,y_2^h,\frac{y_3^h-y_1^h}{2h})\hfill \\ \multicolumn{1}{c}{i=3:\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\frac{y_4^h-2y_3^h+y_2^h}{h^2}}& =\hfill & f(x_3,y_3^h,\frac{y_4^h-y_2^h}{2h})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{with}{y}_4^h=0\hfill \end{array}}$

By directly using the known values at the boundaries there remains a system of equations for ${\overrightarrow{y}}^h=\left(\begin{array}{c}\hfill y_1^h\hfill \\ \hfill :\hfill \\ \hfill y_3^h\hfill \end{array}\right)$. As usual in finite difference methods, the truncation error of this method is defined as the residual one obtains by inserting the true solution $y(x)$ of the boundary value problem into the discretized system, for example in the case of discretizing $y\text{'}(a),\hspace{0.5em}y\text{'}(b)$ by the onesided difference quotient as

${\displaystyle \tau (y;h)\hspace{0.5em}=\hspace{0.5em}\left(\begin{array}{c}\hfill g_1(y(x_0),\frac{y(x_0+h)-y(x_0)}{h},y(x_N),\frac{y(x_N)-y(x_{N-1})}{h})\hfill \\ \hfill \frac{y(x_2)-2y(x_1)+y(x_0)}{h^2}\hspace{0.5em}-\hspace{0.5em}f(x_1,y(x_1),\frac{y(x_2)-y(x_0)}{2h})\hfill \\ \hfill \dots \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\dots \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\dots \hfill \\ \hfill \frac{y(x_N)-2y(x_{N-1})+y(x_{N-2})}{h^2}\hspace{0.5em}-\hspace{0.5em}f(x_{N-1},y(x_{N-1}),\frac{y(x_N)-y(x_{N-1})}{2h})\hfill \\ \hfill g_2(y(x_0),\frac{y(x_0+h)-y(x_0)}{h},y(x_N),\frac{y(x_N)-y(x_{N-1})}{h})\hfill \end{array}\right)}$

which is of order $O(h^2)$ in the middle part and $O(h)$ in the first and the last component. In case of the simple Dirichlet boundary conditions the first and the last component of this vector are zero of course. For this special case considered here there exists an existence and uniqueness theorem, which however requires strong (sufficient) conditions Theorem 1 Assume $f\in C^1([a,b]\times R\times R)$. Moreover assume that there holds

${\displaystyle |\frac{\partial }{\partial s}f(x,r,s)|\le B}$

for all $x\hspace{0.5em}\in \hspace{0.5em}[a,b],\hspace{0.5em}r,s\hspace{0.5em}\in R$ and that

${\displaystyle 0\hspace{0.5em}\le \hspace{0.5em}\frac{\partial }{\partial r}f(x,r,s)}$

for all $x\in [a,b]$, $r,s\in R$.
Then the boundary value problem with homogeneous Dirichlet boundary conditions is uniquely solvable. The discretized system is also uniquely solvable if

${\displaystyle 0<h\le min\{\frac{2}{(B+4{)}^2},b-a\}\hspace{0.5em}.}$

If $f$ is linear in $y$ and $y\text{'}$ then $h\le min\{\frac{2}{B},b-a\}$ suffices for this. The discretized solution ${\overrightarrow{y}}^h$ satisfies

${\displaystyle \oplus \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}|y(x_i)-y_i^h|\le C·h^2}$

for some constant $C$ (i.e. the method has convergence order 2).

If in addition $f\in C^{2m+2}([a,b]\times R\times R)$, then there holds more precisely with some smooth functions $e_i(x)$:

${\displaystyle \oplus \oplus \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}{y}_i^h=y(x_i)+e_1(x_i)h^2+\dots +e_m(x_i)h^{2m}+O(h^{2m+2})\hspace{0.5em}.}$

The last assertion in this theorem says that the global error

${\displaystyle y(x)-y^h(x)}$

behaves itself like a smooth function of $x$ ( $h$ fixed). Now imagine to compute such a discretized solution on grids with size $h,h/2,h/4,\dots$ These grids are coherent and on the crudest grid you win several approximate values of increasing precision. And obviously

${\displaystyle {\tilde{y}}_i^h\hspace{0.5em}\stackrel{\mathrm{def}}{=}\hspace{0.5em}\frac{4y_i^{h/2}-y_i^h}{3}\hspace{0.5em}=\hspace{0.5em}y(x_i)-\frac{1}{4}{e}_2(x_i)h^4+O(h^6)}$

and finally

${\displaystyle \frac{16{\tilde{y}}_i^{h/2}-{\tilde{y}}_i^h}{15}\hspace{0.5em}=\hspace{0.5em}y(x_i)+O(h^6)\hspace{0.5em}.}$

Obviously this technique is the same as used in Romberg's method of quadrature and can be continued to get even higher order of convergence. This is a special case of Richardson's extrapolation to the limit. We use it here for the two stages shown above. Next we have to discuss the case of general boundary conditions. The simplest idea is to discretize here the derivatives by the ordinary difference quotient of first order, getting the two additional equations necessary to fix (locally) the solution.

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{g_1(y_0^h,\frac{y_1^h-y_0^h}{h},y_N^h,\frac{y_N^h-y_{N-1}^h}{h})\hspace{0.5em}}& =\hfill & \hspace{0.5em}0\hfill \\ \multicolumn{1}{c}{g_2(y_0^h,\frac{y_1^h-y_0^h}{h},y_N^h,\frac{y_N^h-y_{N-1}^h}{h})\hspace{0.5em}}& =\hfill & \hspace{0.5em}0\hfill \end{array}}$

This destroys the convergence of second order, leaving a method of order one only. Also, the error expansion in powers of $h$ now contains odd powers, requiring a different and less efficient Richardson extrapolation, since

${\displaystyle \frac{y(x_1)-y(x_0)}{h}\hspace{0.5em}=\hspace{0.5em}y\text{'}(x_0)+\frac{1}{2}hy"(\xi )}$

and similar at the right boundary. There are some tricks to overcome this effect, but they require that the differential equation maintains validity also a bit outside $[a,b]$. That means that the solution can be continued over $[a-\epsilon ,b+\epsilon ]$ for some $\epsilon \hspace{0.5em}>\hspace{0.5em}0$. In this case one introduces for example fictitious nodes $x_{-1}=a-h$ and $x_{N+1}=b+h$ and replaces

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{y\text{'}(a)}& \text{by}\hfill & \frac{y_1^h-y_{-1}^h}{2h}\hfill \\ \multicolumn{1}{c}{y\text{'}(b)}& \text{by}\hfill & \frac{y_{N+1}^h-y_{N-1}^h}{2h}\hfill \end{array}}$

maintaining the development of the global error in even powers of $h$ . Since one has now two unknowns in addition, one uses the discretized differential equation also for $x=a$ and $x=b$ in order to get two additonal equations for them.

3  Solution of the (non)linear system for ${\overrightarrow{y}}^h$

We consider the simple case

${\displaystyle y(a)\hspace{0.5em}=\hspace{0.5em}y(b)\hspace{0.5em}=\hspace{0.5em}0}$

first. Of course in this case

${\displaystyle y_0^h\hspace{0.5em}=\hspace{0.5em}{y}_N^h\hspace{0.5em}=\hspace{0.5em}0}$

too and it remains a system for the unknown values at the $N-1$ interior nodes. The Jacobian of this system then has the tridiagonal form

${\displaystyle \begin{array}{ccc}\multicolumn{1}{c}{}& \hfill & \frac{1}{h^2}(0,\dots ,0,1,-2,1,0,\dots ,0)-(0,\dots ,0,(\frac{\partial }{\partial s}f)(-\frac{1}{2h}),(\frac{\partial }{\partial r}f)·1,(\frac{\partial }{\partial s}f)(\frac{1}{2h}),0,\dots ,0)\hfill \\ \multicolumn{1}{c}{}& =\hfill & \frac{1}{h^2}(0,\dots ,0,1+\frac{h}{2}·(\frac{\partial }{\partial s}f),-2-h^2·(\frac{\partial }{\partial r}f),1-\frac{h}{2}·(\frac{\partial }{\partial s}f),0,\dots ,0)\hspace{0.5em}.\hfill \end{array}}$

The conditions of the above theorem then result in the fact that this Jacobian is always invertible and has a condition number of the order of magnitude of $h^{-2}$. Since the influences of the rounding errors in the computed ${\overrightarrow{y}}^h$ will be in the order of machine precision times condition number and the global discretization error is in the order ${\mathrm{Ch}}^2$ with $C=max\{|e_1(x):\hspace{0.5em}x\hspace{0.5em}\in \hspace{0.5em}[a,b]\}$ it makes little sense to reduce $h$ below the fourth root of the machine precision. Under the strong assumptions of the theorem the solution of the discretized system is globally unique (Hadamards theorem applies) and the damped Newton method will converge from every starting point, but as usual the quadratic convergence will become apparent for good initial guesses only. And one must be aware that the restrictive assumptions of the theorem will usually not apply. At best they will be satisfied in a strip around the true solution and anything will work if the initial guess is in inside this strip. If boundary conditions must also be discretized, then two rows and columns will be added on top and bottom of this matrix. It will maintain the tridiagonal structure if the first order approximation of the derivatives at the boundary are used, otherwise a quindiagonal structure is added there. Anyway the solution of the linearized system within Newton's method will be simple and efficient.

---

File translated from T_EX by T_TM Unregistered, version 4.03.
On 20 Jun 2016, 18:29.