BVP stability¶
AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.
These notebooks are all available on Github.
Compute the inverse of tridiagonal matrices (augmented by boundary conditions) coming from the BVP $u''(x) = f(x)$.
%matplotlib inline
from pylab import *
from scipy.sparse import diags
# import a module with a new local name for brevity:
import scipy.sparse.linalg as sp_linalg
Suppress warnings coming from adding new nonzeros in a crc matrix:
import logging
logging.captureWarnings(True)
Set up the matrix¶
We will investigate the matrix given in equation (2.43) of the textbook for Dirichlet boundary conditions. This uses the formulation where the standard tridiagonal matrix is augmented by two additional rows for $U_0$ and $U_{m+1}$ that correspond to the equations $U_0 = \alpha$ and $U_1 = \beta$.
We will also consider the matrix given at the top of page 32, which is similar but implements a Neumann boundary condition with a second-order accurate one-sided approximation to $u''(x_0)$.
NOTE: The equation at the bottom of page 31 should have $-\sigma$ on the right hand side, and the matrix at the top of page 32 is incorrect. The first row should correspond to the corrected equation from p. 31. The correct version is:
\begin{align*} \frac{1}{h^2} \left[ \begin{array}{ccccccccccccccc} -3h/2 & 2h & -h/2\\ 1&-2&1\\ &1&-2&1\\ &&&\ddots\\ &&&&1&-2&1\\ &&&&&0&h^2 \end{array} \right] ~ \left[ \begin{array}{ccccccccccccccc} U_0 \\ U_1 \\ U_2 \\ \vdots \\ U_m \\ U_{m+1} \end{array} \right] = \left[ \begin{array}{ccccccccccccccc} \sigma \\ f(x_1) \\ f(x_2) \\ \vdots \\ f(x_m) \\ \beta \end{array} \right] \end{align*}
Note that the first equation in this system approximates $u'(x_0) = \sigma$ with a second-order one-sided difference, and the last equation is simply the Dirichlet BC $u(x_{m+1}) = \beta$.
The function below creates such a matrix, and also prints out some information about it, including the norm of the inverse. The bc parameter controls whether the left boundary is Dirchlet or Neumann. In the Dirichlet case the first row is simpler, with only one nonzero.
ax = 0.
bx = 1.
def test_A(m, bc='dirichlet'):
h = (bx-ax) / (m+1)
em = ones(m+2)
em1 = ones(m+1)
A = diags([em1, -2*em, em1], [-1, 0, 1], format='csc')
# fix the first row:
if bc=='dirichlet':
A[0,0] = h**2
A[0,1] = 0.
elif bc=='neumann':
A[0,0] = -3*h/2.
A[0,1] = 2*h
A[0,2] = -h/2. # adding a new nonzero
else:
raise ValueError('Unrecognized bc: %s' % bc)
# fix the last row:
A[m+1,m] = 0.
A[m+1,m+1] = h**2
A = A / h**2
print('m = ', m)
print('A has type %s, of shape %s' % (type(A), A.shape))
Ainv = sp_linalg.inv(A)
normAinv = sp_linalg.norm(Ainv,inf)
print('Infinity norm of Ainv = %g' % normAinv)
return A
Dirichlet boundary conditions¶
Here's what the matrix looks like for a small value of $m$:
A = test_A(5, 'dirichlet')
print(A.toarray())
Note that the max-norm of $A^{-1}$ is 1.125. For stability we hope this is uniformly bounded as we increase $m$ (and decrease $h$). In fact we see it is constant:
A = test_A(99, 'dirichlet')
A = test_A(199, 'dirichlet')
Plot the columns of $A^{-1}$¶
Rather than printing out $A^{-1}$, it is more illuminating to plot the values in each column vs. the row index.
From the discussion of Section 2.11 you should know what values to expect for the interior columns of $A^{-1}$ (see Figure 2.1 in the book). The first and last column are plotted separately below since they are scaled differently. Think about what these columns represent in terms of the way we have included the boundary conditions into the matrix formulation.
m = 5
x = linspace(ax,bx,m+2)
A = test_A(m, bc='dirichlet')
Ainv = sp_linalg.inv(A).toarray()
figure(figsize=(12,5))
subplot(1,2,1)
for j in range(1,m+1):
plot(Ainv[:,j], 'o-', label='column %i' % j)
legend()
xlabel('row index')
ylabel('Ainv value')
subplot(1,2,2)
plot(Ainv[:,0], 'o-', label='column 0')
plot(Ainv[:,m+1], 'o-', label='column %i' % (m+1))
legend()
xlabel('row index')
Neumann boundary conditions¶
Repeat these tests with Neumman conditions:
A = test_A(5, 'neumann')
print(A.toarray())
Note that again the max-norm of $A^{-1}$ stays constant as we increase $m$:
A = test_A(99, bc='neumann')
A = test_A(199, bc='neumann')
Plot the columns of $A^{-1}$¶
Think about why these columns have the form they do.
m = 5
x = linspace(ax,bx,m+2)
A = test_A(m, bc='neumann')
Ainv = sp_linalg.inv(A).toarray()
figure(figsize=(12,5))
subplot(1,2,1)
for j in range(1,m+1):
plot(Ainv[:,j], 'o-', label='column %i' % j)
legend()
xlabel('row index')
ylabel('Ainv value')
subplot(1,2,2)
plot(Ainv[:,0], 'o-', label='column 0')
plot(Ainv[:,m+1], 'o-', label='column %i' % (m+1))
legend()
xlabel('row index')
