Boundary Conditions

Can have tremendous effect on solution. Choice dictated by physics and well-posedness.

Numerical approximations to boundary condition's must be done with great care since they are potentially capable of making an otherwise stable scheme unstable.

There's usually no difficulty implementing Dirichlet type boundary conditions: take boundary conditions at $x=0,$ for simplicity:

$$U(t,0)=b(t)$$$$\mbox { approximated as }u^n_0=b^n$$

No problem with periodic boundary condition's: assume other boundary at $x=L>0:$

$$U(t,0)=U(t,L)$$$$\mbox { approximated as }u^n_0=u^n_{m+1}$$

if there are $M+1$ gridpoints.

What about Neumann? Take $x=0$ and $x=1$ as boundaries:

Consider first $$\frac{\partial U}{d x}\vert _{{x=0}}=a_0$$ and $\frac{\partial U}{\partial x}\vert _{x+1}=a_1$ A first-order approximation would be $$\frac{u^n_1-u^n_0}{h}=a_0$$ and $\frac{u^n_M-U^n_{m-1}}{h}=a_1$.

Using ghost values $u^n_{-1}$ and $u^n_{M+1}$, a second-order approximation

$$\frac{u^n_1-u^n_{-1}}{2h}=a_0$$$$\qquad \mbox {and} \qquad \frac{u^{n}_{M+1}-u_{M-1}}{2h}=a_1$$

One can also use a $2^{nd}$-order approximation, the 1-sided

$$\frac{-3u^n_0+4u^n_1-u^n_2}{2h}=a_0$$

and a similiar expression can be found for the other boundary.

Remarks

• The choice of boundary conditions is dictated by the physics and mathematical well-posedness. However, the discrete version depends on numerical stability considerations. ALWAYS CHECK FOR STABILITY AT ENDPOINTS.
• This is not a theorem, but in general you want to use the same order of finite difference approximation to the boundary conditions as you did for the interior points. What happens if you go with lower order? You defeat the purpose of using a high order scheme. What happens when you go higher order? In general, you get an ill-poosed linear algebraic problem.
• Using one-sided derivatives is ok, but these are prone to generate ill-posed linear algebraic problems if you use mix them. This will be apparent when you do your stability study at the end points.

Note on Advection-Diffusion Equation (UPWINDING IS IMPORTANT)

We're not considering this important equation, other than by an example that illustrates one of the important facets that makes this equation tricky to approximate. We'll take as example the forward- time/central-space scheme.

There are a couple of things we're going to learn: upwinding, and also mixing explicit and implicit methods when you have you're solving a problem of mixed type.

We'll work by example:

$$\frac{u^{n+1}_m-u^n_m}{k}+a\frac{u^n_{m+1}-u^2_{m-1}} {2h}=b\frac{u^n_{m+1}-2u^n_{m}+u^n_{m-1}}{h^2}$$ (140)

as an approximation to

$$u_t+au_x=bu_{xx}, \overbrace {a>0}^{\mbox{for concreteness}} \mbox {and} \overbrace {b>0}^{\mbox {must be}}$$ (141)

not because its a good algorithm, but because it illustrates the point to be made very simply:

Scheme (141) is $(1,2)$ order accurate and $2^{nd}$ order accurate overall because of stability requirement

$$b\mu \le\frac12 {\mbox { (check this!)}}$$

$$\mbox {where }\mu = \frac{h}{k^2}$$

Note that this restriction is rather severe on the time step, since getting good resolution in space by making $h$ small might mean really tiny $k$, which means long computing times.

Now,

$$\begin{array}{ll} \mbox {let }\alpha\equiv \displaystyle \fra... ...r,'' and let $\lambda=\displaystyle \frac{k}{h}$}}. \end{array}$$

The $\alpha$ is a ratio of the importance of inertial (wave-like or signal-propagating effects to diffusion effects). So that if $\alpha<1$ we have a diffusion dominated problem, etc.

Solving for $u^{n+1}_{m}$ and taking the absolute value of both sides yields

$$\vert u^{n+1}_{m}\vert\le (1-2b\mu)\vert u^n_m\vert+b\mu(1-\alpha)\vert u^{n}_{m+1}\vert+b\mu (1+\alpha)\vert u^{n}_{m-1}\vert$$

Since we require that

$$\max_m\vert u_m^{n+1}\vert\le \max_m\vert u^n_m\vert$$ (142)

for parabolic equation approximations, by setting $\alpha\le 1$ we satisfy (143).

Hence we need to satisfy two conditions

$$\left\{ \begin{array}{l} b\mu\le\frac12\\ \alpha\le 1\lef... ...sfied if }} h\le \displaystyle \frac{2b}{a} \end{array}\right.$$

the second is a restriction on mesh spacing and could be very restrictive.

If $b\mu\le\displaystyle \frac12$ is the stability condition, what does $\alpha\le 1$ do? It guarantees that the scheme will behave qualitatively like a parabolic equation approximation. What you will see if you set $\alpha >1$ is the appearance of spurious oscillations $\ldots$ they usually do not grow excessively and they result from inadequate resolution.

One way to avoid the restriction on mesh-spacing is to use ``upwind differencing.'' The scheme then is

$$\frac{u^{n+1}_m-u^n_m}{k}+a\frac{u^n_m-u^n_{m-1}}{h}=b\frac{1}{h^2}\left( u^b_{m+1}-2u^n_m+u^n_{m-1}\right)$$ (143)

or $u^{n+1}_{m}=[1-2b\mu(1+\alpha)]u^n_m+b\mu u^n_{m+1}+b\mu(1+2\alpha)u^n_{m-1}$

If $1-2b\mu(1+\alpha)\rightarrow 0\Rightarrow$ scheme satisfies (143). This is satisfied if

$$2b\mu+a\lambda\le 1$$

Oscillations are eliminated but now we have 1-order accuracy in space. When $b$ small and a large (typical) this condition is much less restrictive than $h\le 2b/a$. Note, however that (144) can be written as

$$\frac{u^{n+1}_m-u^n_m}{k}+a\frac{u^n_{m+1-u^n_{m-1}}}{2h}= \under... ...t(b+\frac{ah}{2}\right)} \Big(\frac{u^n_{m+1}-2u^n_m+u^n_{m-1}}{h^2}\Big) {h^2}$$

$$\begin{array}{c} \mbox {\lq\lq artificial viscosity''}\\ \mbox... ...ed to make solutions}\\ \mbox {non- oscillatory.} \end{array}$$

Remark

• Note that upwinding direction was given by $a$. If $a<0$, the derivative would involve $m+1$ and $m$ values of $u$ instead.
• What if $a$ has different signs for different parts of the domain? Use a switching in your code to check for upwinding direction and then change the derivative to the relative upwind direction.
• One way to avoid oscillations due to constraints on the Peclet number and at the same time get larger time steps for asymptotic stability is to compute the problem by using an explicit upwind in the advective term $U_x$, which will generate a mild restriction on mesh spacing $h$, and then go implicit in the diffusive term $U_xx$, such as Crank Nicholson. The resulting code will be quite robust.