Finite Difference Technique

LINEAR CASE

$$\mbox {Take}\left\{\begin{array}{ll} Y''=p(x)Y'+q(x)Y+r(x) & x_0\le x\le b\\ Y(x_0)=\alpha\\ Y(b)=\beta \end{array}\right.$$

RECIPE

(Equally spaced grid case)

$$\begin{array}{ll} \Rightarrow [x_0, b] & \mbox{ and divide in... ...,1, \cdots N+1\\ h=\displaystyle \frac{(b-x_0)}{N} \end{array}$$

Approximate $Y''$ and $Y'$ by difference quotients

Approximate $Y_n$ by $y_n$

Generate an $N\times N$ matrix problem for the unknowns $y_n$

The boundary data is at $x_0$, at which point $y_0=\alpha$, and at $x_{N+1}$, at which point $y_{N+1} = \beta$.

Example In this instance we'll use low order center-differenced approximations to the derivatives. Assume that $Y\in C^4[x_{i+0}x_{i+1}]$. To get the finite difference expressions to the derivatives, expand

$$Y(x_{i}\pm h)=Y(x_{i\pm 1})=Y(x_i)\pm hY'(x_i)+\frac{h^2}{2}Y''(x_i)\pm$$

$$\frac{h^3}{6}Y'''(x_i)+\frac{h^4}{24}Y''''(\xi_i)$$ (80)

for some $\xi_i^{\pm}$ in $(x_i,x_{i\pm 1})$

Add $Y(x_{i}+h)$ and $Y(x_{i}-h)$ expressions and get

$$Y''(x_i)=\frac{1}{h^2}\left[Y(x_{i+1})-2Y(x_i)+Y(x_{i-1})\right]-\frac{h^2}{24}\left[Y''''(\xi^+_i)+Y''''(\xi^{-}_i \right]$$

Use the Intermediate Value Theorem and (80) to deduce

$$Y''(x_i)=\frac{1}{h^2}\left[Y(x_{i+1}) -2Y(x_i)+Y(x_{i-1})\right]-\frac{h^2} {24}Y''''(\xi_i)$$

with

$$\xi_i\in (x_{i-1}, x_{i+1})$$

So $Y''(x_i)\approx\displaystyle \frac{1}{h^2}\left[Y_{i+1}-2Y_i+Y_{i-1}\right]$. This is a ``centered-difference'' aproximation to $Y''(x_i)$ and is an ${\cal O}(h^2)$

Exercise Show that

$Y'(x_i)\approx \displaystyle \frac{1}{2h}\left[Y_{i+1}-Y_{i-1}\right]$ is $O(h^2)$ approximation to $Y'$ at $x=x_i$. The truncation error is $$\frac{h^2}{6}Y''(\eta_i)$$

$\Box$

So, now we project the equation onto our grid, with $y_i\equiv y(x_i)$, we get

$$\left\{\begin{array}{l} \displaystyle \frac{y_{i+1}-2y_{i}+y_{i-1... ...i+r(x_i)\quad i=1, 2\cdots N y_0=\alpha\quad y_{N+1}=\beta \end{array}\right.$$ (81)

The truncation error is $-\displaystyle \frac{h^2}{12}\left[2p(x_i)y'''(\eta_i)- y''''(\xi_i)\right]={\cal O}(h^2).$ So (81) has the form of the linear system

$$A{\bf y}={\bf b}$$ (82)

with

$$A=\left[\begin{array}{ccccc} 2+h^2q_1 &-1+\displaystyle \frac{h}... ...-1}\\ 0 &0 &0 & -1-\displaystyle \frac{h}{2}p_N &2+h^2q_N \end{array}\right]$$

and

$${\bf y}=\left[\begin{array}{l} y_i\\ :\\ y_N \end{array}\right... ...\ -h^2r_{N-1}\\ -h^2r_N+\left(1-\frac{h}{2}pN\right)\beta \end{array}\right]$$

Theorem

Suppose $p$, $q$, $r$ are continuous on $[x_o, b]$. If $q(x)\ge 0$ on $[x_0,b]\Rightarrow$ (82) has a unique solution provided $h<2/L$ where $L=\displaystyle \max_{x_0\le x\le b}\vert p(x)\vert$. Proof (exercise). Look at conditions for the solution of the associated problem (82). $\Box$ Remark: How could we get higher than ${\cal O}(h^2)$ truncation error? Could use higher order approximation to derivatives, but this leads to more computing (nothing to be scared about these days). However, the more important problem is that it leads to a mathematical issue to resolve: Suppose you use ${\cal O}(h^4)$ approximation to the derivatives: we'll require knowledge of field at $y_{i-2}, y_{i-1}, y_i, y_{i+1} y_{i+2}$. In the interior this is not a problem...our matrix $A$ will simply have a larger bandwidth. However, at the end-points we have some figuring out to do: at $i=1$, we would need $y_{-1}, y_{0}, y_{1}, y_{2}$. The node $y_0$ is given by boundary conditions, $y_1$ and $y_2$ are interior points, so these are ok, but we don't have an equation or condition to determine $y_{1}$. At $i=N$, we have the same sort of problem but there the undertermined node is at $y_{N+2}$. As you might guess, using even a higher order scheme will generate even more undetermined boundary points. So, how do we deal with this problem? If we want to get high order accuracy we need to determine these unknowns.

Strategies:

• Use a right-sided finite difference approximation to derivatives at $i=0$ and a left-sided finite difference approximation to derivatives at $i=N+1$. This seldom works. What you need to do is construct the associated linear algebra problem, (82), and check to see if the system indeed has a solution...what happens in many instances is that you get an underdetermined system.
• Use a lower order discretization and commensurably smaller step sizes at the ends and match the solution there to the interior so that the overall order of the scheme is retained. This is a rather clumsy strategy and sometimes generates very fragile code. It will only work in certain instances. In fact, in some instances, for example, in boundary layer problems, a sudden change in grid size will generate spurious and usually noisy approximations at the end points.
• Use the physics of problem to generate conditions for all undetermined $y_i$'s to the left and the right of $i=1$ and at $i=N$, repectively. This is a very sensible approach, but is not always possible.
• Use Richardson extrapolation (e.g. stick to the second-order approximation and you get higher order but you'll have to compute the problem twice). This is also a very sensible approach and is viable, provided you have precise estimates for the truncation error (not a difficult thing to calculate). Remember, though, that as a rule you don't want to use extrapolation more than a couple of times, since it may become too computationally intensive, as compared to simply using a higher order scheme, or you may encounter rounding-errors. See Richardson Extrapolation.
• Combine the last two strategies.

Unequal Spacing: One other possibility for dealing with the above problem is to use variable-sized meshes. This is not only used to get around the problem of undetermined quantities but is also generally applicable to boundary value problems in which multiple scales in the approximate solution are expected to arise. This is typically used to solve problems which have localized values of $x_i$ for which a lot of changes in the dependent variables are seen, but not much is happening in other areas. This is a topic of current research and is beyond the scope of presentation: in general, one hope to obtain higher resolution in certain places where $y$ varies a lot and go with low resolution in places where $y$ does not change all that much. The overall aim is to reduce the computational expense of the method, as compared to over-resolving the whole domain, but this must be done with care. Perhaps the most research in this area is done by the finite element methods community and if you are interested, you could start by just finding out what finite element methods are. Then look into adaptive mesh refinement schemes and domain decomposition methods.

NONLINEAR CASE:

Take

$$\left\{\begin{array}{l} Y''=f(x,Y,Y') \qquad x_0\le x\le b Y(x_0)=\alpha Y(b)=\beta \end{array}\right.$$ (83)

In principle, it's simple. Same deal as above, but now we get a nonlinear system of algebraic equations. How to solve? Use Newton's method or fixed point iteration. But before launching into either of these two solution techniques, make sure that you study the system you're solving: does the (83) have a solution? Is it unique? Will a fixed point method be applicable? Again, fixed point methods will not require the calculation of an analytical or approximate Jacobian.

The easy ones are those with a unique solution, i.e. when

$$\left\{\begin{array}{ll} \mbox{(i)} & f,f_Y,f_Y' \mbox{ are all c... ...vert f_Y(x,Y,Y')\vert\quad L=\max_D\vert f_{Y'}(x,Y,Y')\vert \end{array}\right.$$ (84)

$\Box$

Recipe: Using ${\cal O}(h^2)$ discretization for $Y$, its derivatives, and the equation coefficients, just as in linear case.

1. Form associated nonlinear algebraic system.
2. Solve system using fixed point methods, or better yet, express the resulting nonlinear system of algebraic equations as

$${\bf S}({\bf y}^{\ell +1})=0$$

with ${\bf y}^{\ell}$= $\left[\begin{array}{l} y_1^{\ell}\\ :\\ .\\ y_{N}^{\ell} \end{array}\right]$, known. The superscript $\ell$ is the iteration counter.

then ${\bf S}({\bf y}^{\ell+1})={\bf S}({\bf y}^{\ell})+J ({\bf y}^{\ell})\delta {\bf y}^{\ell} +{\cal O}(\delta {\bf y^2})=0$ assume ${\cal O} (\delta y^2)$ small, then the approximation

$$J(y^{\ell})\delta y^{\ell} =-{\bf S} (y^{\ell})$$ (85)

where $J({\bf Y}^{\ell})$ is the Jacobian of the nonlinear algebraic problem, evaluated at ${\bf y}^{\ell}$, may be suitable. And $\delta {\bf y}^{\ell}$ is the correction on ${\bf y}^{\ell},$ so that

$${\bf y}^{\ell+1}={\bf y}^{\ell}+\delta {\bf y}^{\ell}$$

Problem (85) is a linear algebraic problem, in fact, it can be shown to be tridiagonal.

Implementation comments: don't forget to put a ``max iterations exceeded'' condition in your code. Also, since scheme is ${\cal O}(h^2)$, use ${\cal O}(h^2)$ stopping criteria for the Newton iteration. You can use fixed point method in the iteration, but convergence is linear and has restrictions on the type of problem for which it is applicable (see Fixed Point iteration. Since a good intial guess is required, make it go through $(x_0,\alpha)$ and $(b,\beta)$ and in addition, satisfy appropriate conditions as per (83) on previous page.

Want to use a higher order truncation scheme? Don't neglect considering a low-order method, coupled with Richardson extrapolation, before getting into something more involved. Otherwise, use a higher-order scheme, but make sure that you don't leave any boundary points undertermined.

$\Box$