BVPs, BCs, nonlinear, nonuniform grids¶
Boundary conditions¶
- Dirichlet conditions are straightforward
- Unknowns are interior points. Points near the boundaries are written in terms of the boundaries.
- Neumann and Robin conditions:
- $y^{\prime} = \alpha$ is given at the boundary (or $y^{\prime} + \beta y = \alpha$ for Robin).
- Interior cells need $y$ on the boundary, not $y^{\prime}$.
Two Approaches¶
Approach 1. Ghost Cell Method Include the boundary point in the list of unknowns (like interior points).
- Discretize the boundary point exactly as if it was an interior point.
- For central differences, this will reference one point past the boundary.
\\| \\| \\| * * * * i-1 \\|i i+1 -1 \\|0 1 2 \\| \\|
- Use the BC $y^{\prime}=\alpha$ as an equation for the new unknown $y_{-1}$. $$y_0^{\prime} = \alpha \rightarrow \frac{y_1-y_{-1}}{2\Delta x} = \alpha \rightarrow$$ $$y_{-1}=y_{1} - 2\Delta x\alpha.$$
- Substitute this equation in for $y_{i-1} = y_{-1}$ appearing in the finite difference equation (FDE) at point $i=0$. The first two equations below are the FDE at $i=0$. The last equation has the substitution for $y_{-1}$: $$l_iy_{i-1} + a_iy_i + u_{i}y_{i+1} = F_i$$ $$\mbox{or,}$$ $$l_0y_{-1} + a_0y_0 + u_0y_{1} = F_0, $$ $$\rightarrow a_0y_0 + (l_0+u_0)y_1 = F_0 + 2\Delta x\alpha l_0. $$
- That is, the unknown $y_{-1}$ is in terms of $y_1$, so when we substitute it into the FDE at $i=0$, the coefficient of $y_1$ and the RHS are modified in the $i=0$ equation.
- This is called the Ghost Cell Method since a false or ghost point arises, which is then handled with the boundary condition equation.
Advantages
- Solves directly for the unknown boundary value.
- Uses a uniform stencil (central difference everywhere, even at the boundaries.
Disadvantage
- Higher order $y^{\prime\prime\prime}$ or higher order FDA can lead to multiple outside points, which can be awkward.
- Remedy this using one-sided differences.
Method 2 Don't include the boundary point in the list of unkowns.
```
\\|
\\|
\\|
* * * *
\\|i-1 i i+1
\\|-1 0 1 2
\\|
```
- The first unknown point is $i=0$.
- The FDE at this point is in terms of point $i-1=-1$
- Use the Boundary condition at $i-1$ with a one-sided difference to get another equation to write $y_{-1}$ in terms of the interior points we are solving for ($y_0$, $y_{1}$, etc.)
- At the left side for the $i=0$ FDE (which has a $y_{-1}$ term):
- Write a one sided difference for $y^{\prime}_{-1}$. $$y^{\prime}_{-1} = \alpha = \frac{-\frac{3}{2}y_{-1} + 2y_0 - \frac{1}{2}y_1}{\Delta x}.$$
- Solve for $y_{-1}$: $$ y_{-1} = \frac{\alpha\Delta x - 2y_0 + \frac{1}{2}y_{1}}{-3/2}.$$
- The FDE at point $i=0$ is $$l_0y_{-1} + a_0y_0 + u_0y_1 = F_0.$$
- Insert the above BC equation for $y_{-1}$ (blue) into this FDE to get the final FDE at the $i=0$ point. $$ (a_0 + \frac{4}{3}l_0)y_0 + (u_0 - \frac{1}{3}l_0)y_1 = F_0 + \frac{2}{3}\alpha l_0\Delta x.$$
- A similar procedure is done at the right side of the domain.
- Advantage
- Extends easily to higher order
- On the homework, if we exclude the BC point, then we don't have to divide by $r=0$ at the cylindrical centerline.
- Disadvantage
- Two different stencils are needed. One for the fully interior points $i=1,\,2,\ldots$, and one for the point $i=0$ next to the boundary. (Similarly for the right side of the domain.)
- The boundary point remains unknown.
- Once the solution to the interior points $i=0,\,1,\ldots$ is found, we can find $y_{-1}$, (the boundary point), using the above blue equation.
Nonlinear relaxation methods¶
- The usual approach is to iterate
- Use a linearized form of the equation.
- Guess a solution.
- Iterate to improve it.
- Apply the FDA to $y^{\prime\prime}$, $y^{\prime}$ as before.
- But consier a term like $y^\prime y$
- This results in nonlinear terms: $$y^{\prime}y \rightarrow \left(\frac{y_{i+1}-y_{i-1}}{2\Delta x}\right)\cdot y_i.$$
- We can linearize these terms by splitting the product and lagging part of it. For example $$y_i^2\rightarrow y_i\cdot y_i \rightarrow y_i^{new}\cdot y_i^{old},$$ where $y_i^{old}$ is the value from the previous iteration, which is known.
- For $$y^{\prime\prime} + P(x,y)y^{\prime} + Q(x,y)y = F(x),$$ use $$y^{\prime\prime} + P(x,y^{old})y^{\prime} + Q(x,y^{old})y = F(x),$$
- Example
$$\nabla\cdot\vec{q} = F \rightarrow -\nabla\cdot(k\nabla T) = F.$$
$$-\frac{d}{dx}\left(k\frac{dT}{dx}\right) = F.$$
$$\frac{dT}{dx}\frac{dk}{dx} + k\frac{d^2T}{dx^2} = -F.$$
- $k = k(T)$ so the above equation is nonlinear. So lag $T$ when evaluating $k$. Use $k(T^{old})$.
- This is now linear, but we have to iterate.
Question What are other approaches that could (or should?) be used?
- Newton's method?
- What about a Taylor Series linearization instead of the linearization shown?
- Consider $yy$
- above linearization: $yy\approx y_0y$
- Taylor series: $yy\approx y_0^2 + 2y_0(y-y_0)$
- Plot these for some arbitrary $y_0$
- Consider $yy$
In [14]:
import numpy as np
import matplotlib.pyplot as plt
y = np.linspace(0,10,100)
y0 = 3
y2 = y**2
y0y = y0 * y
yta = y0**2 + 2*y0*(y-y0)
plt.rc('font', size=14)
plt.plot(y,y2, color='grey', lw=5)
plt.plot(y,yta, 'k-')
plt.plot(y,y0y, 'b--')
plt.xlabel("y")
plt.ylabel("yy")
plt.legend(['exact', 'Taylor', 'simple'], frameon=False);
Nonuniform grids¶
Method 1¶
- Set an arbitrary grid with arbitrary $\Delta x$ spacing between points.
* * *
i-1 i i+1
- Let $\Delta x_{i-1} = x_i - x_{i-1}$, and $\Delta x_{i} = x_{i+1}-x_i$.
- Then a central difference approximation is
$$f^{\prime}(x_i) = f_i^{\prime} \approx \frac{f_{i+1}-f_{i-1}}{\Delta x_{i-1} + \Delta x_i}.$$
- A Taylor Series gives second order when $\Delta x_{i-1}=\Delta x_i$, but only first order when $\Delta x_{i-1}\ne \Delta x_i$.
- These are second and first order asymptotically as $\Delta x\rightarrow 0$.
- In practice, accuracy is not severely compromized if $\Delta x_{i-1}$ is not too different from $\Delta x_i$.
- A Taylor Series gives second order when $\Delta x_{i-1}=\Delta x_i$, but only first order when $\Delta x_{i-1}\ne \Delta x_i$.
- For a second derivative, we have $$f^{\prime\prime}_i\approx \frac{\frac{f_{i+1}-f_i}{\Delta x_i}-\frac{f_i-f_{i-1}}{\Delta x_{i-1}}}{\frac{\Delta x_{i-1}+\Delta x_i}{2}}.$$
Method 2¶
- Stretch the grid analytically
- Let $x$ be the nonuniform grid and
- Let $\eta$ be a corresponding uniform grid.
- For example $$\eta = \ln(x+1) \rightarrow x = e^{\eta} - 1.$$
- Then
$$ \frac{dy}{dx} = \frac{dy}{d\eta}\frac{d\eta}{dx} = \frac{dy}{d\eta}\left(\frac{1}{x+1}\right).$$
- This allows us to transform all the $dy/dx$ derivatives on the nonuniform grid to $dy/d\eta$ derivatives on a uniform grid times some known factor $1/(x+1)$ at any given point.
- So, if we are evaluating $dy/dx$ at point $i$, we would have $$y^{\prime}_i = \frac{y_{i+1}-y_{i-1}}{\Delta \eta}\left(\frac{1}{x_i+1}\right).$$
- This works great, and is easy to use.
Questions¶
- What would you do to have an arbitrary grid spacing that is not defined by some known function?
- Could you develop a second order approximation to a central derivative on a nonuniform grid?
In [ ]: