"""Demo of multigrid features for :file:`multigrid.rst`."""
from __future__ import division
__all__ = ['Problem', 'Galerkin', 'MultigridBase', 'MultigridPeriodic',
'MultigridNeumann', 'comparing_periodic']
import math
import numpy as np
import scipy.integrate
import scipy as sp
from mmf.objects import StateVars, Container, process_vars
from mmf.objects import Computed, Delegate
[docs]class Problem(StateVars):
r"""Representation of the following linear problem for use with a
one-dimensional multigrid method. We provided the following
problem:
.. math::
\begin{aligned}
-f''(x) &= a\pi^2\left[\cos(\pi x)-a\sin^2(\pi x)\right]
e^{a(\cos\pi x -1)}, \\
-f''(x) + a^2\pi^2\sin^2(\pi x)f(x) &=
a\pi^2\cos(\pi x)e^{a(\cos\pi x -1)}, \\
-f''(x) - a\pi^2\cos(\pi x)f(x) &=
-a^2\pi^2\sin^2(\pi x)e^{a(\cos\pi x -1)}, \\
f(x) &= e^{a(\cos\pi x - 1)}.
\end{aligned}
This can also be written as a non-linear equation:
.. math::
-f''(x) = \pi^2\left[\ln f + a - a^2\sin^2(\pi x)\right]f
We may also write this as a set of coupled equations in a few
ways, including:
.. math::
\begin{aligned}
- f''(x) &= a\pi^2[g - a(1-g^2)]f, \\
- g''(x) &= \pi^2 g = \pi^2 (a^{-1}\ln f + 1)
\end{aligned}
This gives rise to the following systems for the Jacobian:
.. math::
G_3(f, g) = \begin{pmatrix}
- f''(x) - a\pi^2[g - a(1-g^2)]f, \\
- g''(x) - \pi^2 g
\end{pmatrix}, \qquad
J_3 = -\diff[2]{}{x} + \begin{pmatrix}
- a\pi^2[g - a(1-g^2)] & - a\pi^2[1 + 2g]f\\
0 & - \pi^2
\end{pmatrix}.
.. math::
G_4(f, g) = \begin{pmatrix}
- f''(x) - a\pi^2[g - a(1-g^2)]f, \\
- g''(x) - \pi^2 (a^{-1}\ln f + 1)
\end{pmatrix}, \qquad
J_4 = -\diff[2]{}{x} + \begin{pmatrix}
- a\pi^2[g - a(1-g^2)] & - a\pi^2[1 + 2g]f\\
- \pi^2 a^{-1}/f & 0
\end{pmatrix}.
"""
_state_vars = [
('a', 5,
'The larger this is, the more sharply peaked the solution.'),
('periodic', True,
'''If `True`, then periodic boundary conditions are
represented with $L=2$ ranging from $x\in [0, 2]$,
otherwise Neumann boundary conditions are represented by
with $L=1$ and $x\in [0, 1]$.'''),
('prob', 0,
'''Value in `[0, 1, 2, 3, 4]` specifies which of the three
problems is represented. These differ by which terms are
included in the equation and which terms are included in
the rhs:
`0` : represents the simple Poisson problem
`1` : includes the $\sin^2$ piece which should make the problem
somewhat easier (being positive).
`2` : includes the $\cos$ piece (which is challenging
since it is negative)
`3` : represents the Newton step for problem $G_3$.
`4` : represents the Newton step for problem $G_4$.
'''),
('L', Computed, 'Length'),
]
process_vars()
def __init__(self, *v, **kw):
if self.periodic:
self.L = 2
else:
self.L = 1
def D(self, f, x):
r"""Return $D(f)$, the diagonal portion of the problem (not including
the operator $-T = -\nabla^2$) acting on `f`."""
a = self.a
pi = np.pi
c = np.cos(pi*x)
s = np.sin(pi*x)
if 0 == self.prob:
# Simple problem
return 0*x
elif 1 == self.prob:
return f*pi**2*(a*s)**2
elif 2 == self.prob:
return -f*pi**2*(a*c)
else:
raise NotImplementedError
def DG(self, fg, x):
r"""Return the objective function for Newton's method (without
the $-T$ portion)."""
N = len(x)
a = self.a
pi = np.pi
G = 0*fg
fg[:N] = -a*pi**2*((g - a*(1-g**2))*f)
if 3 == self.prob:
G[N:] = -pi**2*g
elif 4 == self.prob:
G[N] = -pi**2*(np.log(f)/a + 1)
else:
raise NotImplementedError
return G
def DJ(self, dfg, f, g, x):
r"""Return $D_J(f, g)(d_f, d_g)$, the block diagonal portion of the
Jacobian (not including the operator $-T = -\nabla^2$) at
`(f, g)` acting on the step `(df, dg)`."""
N = len(x)
df = dfg[:N]
dg = dfg[N:]
a = self.a
pi = np.pi
jfg = 0*dfg
jfg[:N] = -a*pi**2*((g - a*(1-g**2))*df + (1 + 2*g)*f*dg)
if 3 == self.prob:
jfg[N:] = -pi**2*dg
elif 4 == self.prob:
jfg[N:] = -pi**2*df/a/f
else:
raise NotImplementedError
jfg *= -1
return jfg
def b(self, x):
r"""Return right hand side."""
f = np.exp(self.a*(np.cos(np.pi*x) - 1))
a = self.a
pi = np.pi
c = np.cos(pi*x)
s = np.sin(pi*x)
if 0 == self.prob:
return f*pi**2*(a*c - (a*s)**2)
elif 1 == self.prob:
return f*pi**2*(a*c)
elif 2 == self.prob:
return f*pi**2*(-(a*s)**2)
else:
raise NotImplementedError
def f(self, x, d=0):
r"""Return exact solution to continuum problem and its
derivatives (for checking discretization errors)."""
f = np.exp(self.a*(np.cos(np.pi*x) - 1))
a = self.a
pi = np.pi
c = np.cos(pi*x)
s = np.sin(pi*x)
if 0 == d:
return f
elif 1 == d:
return f*pi*(-a*s)
elif 2 == d:
return f*pi**2*(-a*c + (a*s)**2)
elif 3 == d:
return f*pi**3*(3*a**2*c*s + (-a*s)**3 + a*s)
elif 4 == d:
return f*pi**4*(-6*a**3*c*s**2 + (a*s)**4
- 4*(a*s)**2 + 3*(a*c)**2
+ a*c)
else:
raise NotImplementedError
def g(self, x, d=0):
r"""Return exact solution to continuum problem and its
derivatives (for checking discretization errors)."""
a = self.a
pi = np.pi
c = np.cos(pi*x)
s = np.sin(pi*x)
if 0 == d:
return c
elif 1 == d:
return -pi*s
elif 2 == d:
return -pi**2*c
elif 3 == d:
return pi**3*s
elif 4 == d:
return pi**4*c
else:
raise NotImplementedError
[docs]class Galerkin(object):
r"""Provides a Galerkin operator A by interpolating and then
restricting to the top level. This is not efficient
computationally, but can be used to test convergence."""
[docs] def A(self, f):
"Apply the operator Galerkin restriction of the operator A"
n = 0
while len(f) < self.N:
f = self.I(f)
n += 1
f = -self.T(f) + self.D(f, self.x(len(f)))
while n > 0:
f = self.R(f)
n -= 1
return f
[docs]class MultigridBase(StateVars, Galerkin):
r"""Common base for :class:`MultigridNeumann` and
:class:`MultigridPeriodic`."""
_state_vars = [
('n_0', 2, "Minimum size of grid when `n=0`"),
('n', 7, "Number of levels above `n_0`"),
('dx_n', Computed, "Grid spacing of coarsest grid"),
('N', Computed, "Number of grid points"),
('neutralize', False,
r"""In the case of the pure Possion equation, the matrix
`A` is singular and the constant piece of the solution
should be removed. If `True`, then this is done at the
end of :meth:`v_cycle` by calling :meth:`make_neutral`."""),
('smoothing_method', 0,
r"""Smoothing method. One of:
`0` : Weighted Jacobi Iteration
`1` : Red-black Gauss-Seidel
`2` : Lexicographic Gauss-Seidel
"""),
('n_smooth', [2, 2],
r"""`[n_before, n_after]`. Number of times the smoothing
operator is applied before and after the iterative call in
:meth:`v_cycle`."""),
('problem', Delegate(Problem, ['a', 'prob'])),
]
process_vars()
def __init__(self, *v, **kw):
self.N = self.get_N(self.n)
self.dx_n = np.diff(self.x()[0:2])[0]*2**self.n
def get_N(self, n, inverse=False):
r"""Return the number of points in the `n`'th lattice."""
if inverse:
return np.log2((n - 1)/(self.n_0 - 1))
else:
return 2**n*(self.n_0 - 1) + 1
def residue(self, f, b):
r"""Compute the residue $\vect{b} - \mat{A}\vect{f}$
.
"""
res = b - self.A(f)
return res
def v_cycle(self, f, b):
r"""V-cycle.
Parameters
----------
f : array
Solution to be smoothed.
b : array
Right hand side.
"""
m = int(self.get_N(len(f), inverse=True))
if m == 0:
A_ = self.to_mat(self.A, len(f))
f = np.linalg.lstsq(A_, b)[0]
else:
for n in xrange(self.n_smooth[0]):
f = self.smooth(f, b)
err = -self.R(self.residue(f, b))
f -= self.I(self.v_cycle(0*err, err))
for n in xrange(self.n_smooth[1]):
f = self.smooth(f, b)
if self.neutralize:
f = self.make_neutral(f)
return f
def f_cycle(self, f, b):
r"""F-cycle."""
res = self.residue(f, b)
rs = [res]
for i_ in reversed(xrange(self.n)):
rs.insert(0, self.R(rs[0]))
corr = self.v_cycle(0*rs[0], rs[0])
for m in xrange(1, self.n+1):
corr = self.v_cycle(self.I(corr), rs[m])
return f + corr
def fmg(self, f=None):
r"""Full multigrid cycle."""
b_n = self.problem.b(self.x())
if f is not None:
f = 0*b_n
b_n = self.residue(f, b_n)
bs = [b_n]
for i_ in reversed(xrange(self.n)):
bs.insert(0, self.R(bs[0]))
corr = self.v_cycle(0*bs[0], bs[0])
for m in xrange(1, self.n+1):
corr = self.v_cycle(self.I(corr), bs[m])
if f is None:
return corr
else:
return f + corr
def to_mat(self, op, N):
zero = np.zeros(N, dtype=float)
zero[0] = 1
f = op(zero)
zero[0] = 0
mat = np.zeros((len(f), N), dtype=f.dtype)
mat[:, 0] = f
for i_ in xrange(1, N):
zero[i_] = 1
mat[:, i_] = op(zero)
zero[i_] = 0
return mat
def mats(self):
r"""Return a Container of lists of the various matrices."""
xs = []
As = []
Rs = []
Is = []
Ss = []
Ts = []
for m in xrange(self.n+1):
M = self.get_N(m)
xs.append(self.x(M))
if 0 == m:
# May not be a restriction from the smallest lattice
Rs.append(None)
else:
Rs.append(self.to_mat(self.R, M))
Is.append(self.to_mat(self.I, M))
As.append(self.to_mat(self.A, M))
Ts.append(self.to_mat(self.T, M))
Ss.append(self.to_mat(
lambda x:self.smooth(x, b=0*xs[-1]), M))
bs = [self.problem.b(self.x())]
f_ = np.linalg.lstsq(As[-1], bs[0])[0]
if self.neutralize:
f_ = self.make_neutral(f_)
fs = [f_]
for i_ in reversed(xrange(self.n)):
bs.insert(0, self.R(bs[0]))
m = int(self.get_N(len(bs[0]), inverse=True))
f_ = np.linalg.lstsq(As[m], bs[0])[0]
if self.neutralize:
f_ = self.make_neutral(f_)
fs.insert(0, f_)
return Container(x=xs, A=As, R=Rs, I=Is, S=Ss, T=Ts, b=bs, f=fs)
[docs]class MultigridPeriodic(MultigridBase):
r"""Multigrid class for use with Periodic boundary conditions. The
interpolation and restriction operators here are:
The restriction matrix for $n_0=2$ and $n=2\rightarrow 1$ is
.. math::
\mat{R} = \tfrac{1}{2}\mat{I}^{T} = \frac{1}{2}\begin{pmatrix}
1 & \tfrac{1}{2} & 0 & 0 & 0 & 0 & 0 & \tfrac{1}{2} \\
0 & \tfrac{1}{2} & 1 & \tfrac{1}{2} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & \tfrac{1}{2} & 1 & \tfrac{1}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \tfrac{1}{2} & 1 & \tfrac{1}{2}
\end{pmatrix}
The interpolation matrix for $n_0=2$ and $n=1\rightarrow 2$ is
.. math::
\mat{I} = \begin{pmatrix}
1 & 0 & 0 & 0 \\
\tfrac{1}{2} & \tfrac{1}{2} & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & \tfrac{1}{2} & \tfrac{1}{2} & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & \tfrac{1}{2} & \tfrac{1}{2} \\
0 & 0 & 0 & 1\\
\tfrac{1}{2} & 0 & 0 & \tfrac{1}{2} \\
\end{pmatrix}
which preserve the Galerkin condition for operators of the form $T$
(yet to check)
.. math::
\mat{T}_{\text{Periodic}} = \frac{1}{\delta^2}\begin{pmatrix}
-2 & 1 & 0 & 1 \\
1 & -2 & 1 & 0 \\
0 & 1 & -2 & 1 \\
1 & 0 & 1 & -2 \\
\end{pmatrix}.
in the sense that $\mat{R}\mat{T}\mat{I} \simeq \mat{T}/4$ has the
same tridiagonal form. In addition, the operator satisfies the
symmetry and charge-neutrality conditions:
.. math::
\mat{T} = \mat{T}^{T}, \qquad \sum_{i}\mat{T}_{ij} = 0.
One has the following sequence of grid sizes:
.. math::
N = 2^{n}n_0
"""
_state_vars = [
('problem.periodic', True),
]
process_vars()
def __init__(self, *v, **kw):
self.N = self.get_N(self.n)
self.dx_n = np.diff(self.x()[0:2])[0]*2**self.n
def x(self, N=None):
r"""Return abscissa for grid.
The abscissa are:
.. math::
x_{i} = \delta i, \qquad
\delta = \frac{L}{N+1}, \qquad i=0, ..N-1
::
0 x x x x ... x x o
|----|----|----|----|-- ... |----|----|
-d 0 d 2d 3d ... L-d 0 = x
-1 0 1 2 3 ...N-2 N-1 0 = N
R B R B B R
"""
if N is None:
N = self.N
if self.inclusive:
raise(NotImplementedError)
else:
#L = 2 defined in self.problem.L
dx = self.problem.L/(N)
x0 = 0
return np.linspace(x0, self.problem.L - x0, (N+1))[0:-1]
def make_neutral(self, f):
r"""Subtract the constant piece from `f` using the appropriate
integration and subtraction scheme.
If :attr:`exclusive`, then
.. math::
\int_0^L f(x) \d{x} \approx \sum_{n=0}^{N-1} f_n
and the constant constributes as $\int c = cN$.
"""
f_sum = sum(f)
if self.inclusive:
Raise(NotImplementedError)
else:
c = len(f)
return f - f_sum/c
def R(self, f):
r"""Restriction operator:
.. math::
\mat{R} = \tfrac{1}{2}\mat{I}^{T} = \frac{1}{2}\begin{pmatrix}
1 & \tfrac{1}{2} & 0 & 0 & 0 & 0 & 0 & \tfrac{1}{2} \\
0 & \tfrac{1}{2} & 1 & \tfrac{1}{2} & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & \tfrac{1}{2} & 1 & \tfrac{1}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \tfrac{1}{2} & 1 & \tfrac{1}{2}
\end{pmatrix}
"""
res = 1*f[0::2]
res = res+f[1::2]/2.0
res = res+np.append(f[-1], f[1:-1:2])/2.0
res /= 2
return res
def I(self, f):
r"""Interpolation operator:
.. math::
\mat{I} = \begin{pmatrix}
1 & 0 & 0 & 0 \\
\tfrac{1}{2} & \tfrac{1}{2} & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & \tfrac{1}{2} & \tfrac{1}{2} & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & \tfrac{1}{2} & \tfrac{1}{2} \\
0 & 0 & 0 & 1\\
\tfrac{1}{2} & 0 & 0 & \tfrac{1}{2} \\
\end{pmatrix}
"""
N = len(f)
N_new = 2*N
res = np.empty(N_new, dtype=float)
res[0::2] = f
res[1::2] = (f[0::1] + np.append(f[1::1], f[0]))/2.0
return res
def T(self, f, diag=False):
r"""Act with the operator $T(f)=f''$ at the highest level.
.. math::
\mat{T}_{\text{Periodic}} = \frac{1}{\delta^2}\begin{pmatrix}
-2 & 1 & 0 & 1 \\
1 & -2 & 1 & 0 \\
0 & 1 & -2 & 1 \\
1 & 0 & 1 & -2 \\
\end{pmatrix}.
.. note:: In order to preserve the Galerkin properties, we must
ensure that the step size $\delta_n = 2\delta_{n-1}$.
Parameters
----------
f : array
diag : bool
If `True`, then return the diagonals.
"""
m = self.get_N(len(f), inverse=True)
dx = self.dx_n/2**m
if diag:
res = -2*np.ones(len(f), dtype=f.dtype)
else:
f_aug = np.empty(len(f)+2, dtype=f.dtype)
f_aug[1:-1] = f
f_aug[0] = f[-1]
f_aug[-1] = f[0]
res = -2*f + f_aug[2:] + f_aug[:-2]
res /= dx**2
return res
def D(self, f, x):
r"""Return $D(f)$, the diagonal portion of the problem (not including
the operator $-T = -\nabla^2$) acting on `f`.
This simply calls ``self.problem.D``. It is intended that
this only be called on the top-level problem.
"""
return self.problem.D(f, x)
def smooth(self, f, b):
r"""Smoothing operator. The method depends on :attr:`smoothing_method`:
`0` : Weighted Jacobi:
.. math::
\vect{f} \rightarrow \vect{f} - \tfrac{2}{3}\mat{d}^{-1}
(\mat{A}\vect{f} - \vect{b})
`1` : Red-Black Gauss-Seidel:
If :attr:`inclusive` is `False`, then the abscissa are
interior:
.. math::
x_{n} = \frac{\delta}{2} + \delta n, \qquad
\delta = \frac{L}{N}
::
0 x x x x ... x x o
|----|----|----|----|-- ... |----|----|
-d 0 d 2d 3d ... L-d 0 = x
-1 0 1 2 3 ...N-2 N-1 0 = N
R B R B B R
Implemented only for weight = 1
.. math::
\vect{z}^{m+1}_{red} = \mat{d}^{-1}
(\vect{f}^{m} - A_{black} \vect{f}^{m}_{black})\\
\vect{z}^{m+1}_{black} = \mat{d}^{-1}
(\vect{f}^{m} - A_{red} \vect{f}^{m}_{red})\\
\vect{z}^{m+1}_{black} = \mat{d}^{-1}
(\vect{f}^{m} - A_{red} \vect{f}^{m}_{red})\\
\vect{z}^{m+1}_{red} = \mat{d}^{-1}
(\vect{f}^{m} - A_{black} \vect{f}^{m}_{black})
`2` : Lexicographic Gauss-Seidel:
.. math::
\vect{z}^{m+1}_{red} = \mat{d}^{-1}
(\vect{b} - \mat{L}\vect{f}^{m+1} - \mat{U}\vect{f}^{m})\\
\vect{f}^{m+1} = \vect{f}^{m}
+ weight*(\vect{z}^{m+1}-\vect{f}^{m})
\vect{z}^{m+2}_{red} = \mat{d}^{-1}
(\vect{b} - \mat{U}\vect{f}^{m+2} - \mat{L}\vect{f}^{m+1})\\
\vect{f}^{m+2} = \vect{f}^{m+1}
+ weight*(\vect{z}^{m+2}-\vect{f}^{m+1})
Parameters
----------
f : array
Solution to be smoothed.
b : array
Right hand side.
"""
m = self.get_N(len(f), inverse=True)
res = -self.residue(f, b)
if 0 == self.smoothing_method:
# Weighted Jacobi
d = (-self.T(f, diag=True)
+ self.problem.D(0*f + 1, self.x(len(f))))
return f - 2/3*res/d
elif 1 == self.smoothing_method:
#weighted Red-Black Gauss-Seidel
#f[red] = d**(-1)(b[red]-A_black*f[black])
# = d**(-1)(b[red]+T_black*f[black])
# = d**(-1)(b[red]+f[black]/dx**2)
#and similarly for black
weight = 1
if weight != 1:
raise(NotImplementedError)
if ((f.shape[0]==0) |(np.mod(f.shape[0], 2)==1)):
raise(NotImplementedError)
dx = self.dx_n/2**m
d = (-self.T(f, diag=True)
+ self.problem.D(0*f + 1, self.x(len(f))))
f[0::2] = (b[0::2] + f[1::2]/dx**2
+ np.append(f[-1], f[1:-1:2])/dx**2)/d[0::2]
f[1::2] = (b[1::2] + f[0::2]/dx**2
+ np.append(f[2::2], f[0])/dx**2)/d[1::2]
f[1::2] = (b[1::2] + f[0::2]/dx**2
+ np.append(f[2::2], f[0])/dx**2)/d[1::2]
f[0::2] = (b[0::2] + f[1::2]/dx**2
+ np.append(f[-1], f[1:-1:2])/dx**2)/d[0::2]
#we operate twice to maintain symmetry of the operator.
return f
else:
#weighted Lexicographic Gauss-Seidel
weight = 1
if ((f.shape[0]==0)):
raise(NotImplementedError)
d = (-self.T(f, diag=True)
+ self.problem.D(0*f + 1, self.x(len(f))))
dx = self.dx_n/2**m
f[0] = (1-weight)*f[0] + weight*(b[0]+f[1]/dx**2+f[-1]/dx**2)/d[0]
for i in range(1, f.shape[0]-1):
f[i] = (1-weight)*f[i] +\
weight*(b[i]+f[i+1]/dx**2+f[i-1]/dx**2)/d[i]
f[-1] = (1-weight)*f[-1] +\
weight*(b[-1]+f[-2]/dx**2+f[0]/dx**2)/d[-1]
return f
[docs]class MultigridNeumann(MultigridBase):
r"""Multigrid class for use with Neumann boundary conditions. The
interpolation and restriction operators here are:
.. math::
\mat{I} = \begin{pmatrix}
1\\
\tfrac{1}{2} & \tfrac{1}{2}\\
& 1 \\
& \tfrac{1}{2} & \tfrac{1}{2}\\
& & \ddots\\
& & & & \tfrac{1}{2} & \tfrac{1}{2}\\
& & & & & 1
\end{pmatrix} \qquad
\mat{R} = \tfrac{1}{2}\mat{I}^{T} = \frac{1}{2}\begin{pmatrix}
1 & \tfrac{1}{2}\\
& \tfrac{1}{2} & 1 & \tfrac{1}{2}\\
&& \ddots\\
& & & & \tfrac{1}{2} & 1 & \tfrac{1}{2}\\
& & & & & & 1
\end{pmatrix}
which preserve the Galerkin condition for operators of the form
.. math::
\mat{T}_{\text{Neumann}} = \frac{1}{\delta^2}\begin{pmatrix}
-1 & 1 & \\
1 & -2 & 1\\
& \ddots & \ddots & \ddots\\
& & 1 & -2 & 1\\
& & & 1 & -1 \\
\end{pmatrix}.
in the sense that $\mat{R}\mat{T}\mat{I} \simeq \mat{T}/4$ has the
same tridiagonal form. In addition, the operator satisfies the
symmetry and charge-neutrality conditions:
.. math::
\mat{T} = \mat{T}^{T}, \qquad \sum_{i}\mat{T}_{ij} = 0.
One has the following sequence of grid sizes:
.. math::
N = 2^{n}(n_0 - 1) + 1
"""
_state_vars = [
('inclusive', False,
r"""One can use either interpretation of the boundary
conditions. If :attr:`inclusive` is `True`, then the abscissa
include the endpoints:
.. math::
x_{n} = \delta n, \qquad
\delta = \frac{L}{N-1}
::
x x x ... x x
|---|---|-- ... |---|
0 d 2d ... L
Note that in this case, the proper integration rule is the
trapezoidal rule with factors of 1/2 at the endpoints (see
:meth:`make_neutral`):
.. math::
\int_0^L f(x) \d{x} \approx
\frac{f_0}{2} + f_1 + \cdots + f_{N-2} + \frac{f_{N-1}}{2}
If :attr:`inclusive` is `False`, then the abscissa are
interior:
.. math::
x_{n} = \frac{\delta}{2} + \delta n, \qquad
\delta = \frac{L}{N}
::
o x x x ... x x o
|--|----|----|-- ... |----|--|
0 d/2 3d/2 5d/2 ... L
The proper integration rule is now a midpoint trapezoidal
rule without factors of 1/2 (see :meth:`make_neutral`):
.. math::
\int_0^L f(x) \d{x} \approx \sum_{n=0}^{N-1} f_n
"""),
('problem.periodic', False),
]
process_vars()
def __init__(self, *v, **kw):
self.N = self.get_N(self.n)
self.dx_n = np.diff(self.x()[0:2])[0]*2**self.n
def x(self, N=None):
r"""Return abscissa for grid."""
if N is None:
N = self.N
if self.inclusive:
dx = self.problem.L/(N-1)
x0 = 0
else:
dx = self.problem.L/N
x0 = dx/2
return np.linspace(x0, self.problem.L - x0, N)
def make_neutral(self, f):
r"""Subtract the constant piece from `f` using the appropriate
integration and subtraction scheme.
If :attr:`inclusive`, then
.. math::
\int_0^L f(x) \d{x} \approx
\frac{f_0}{2} + f_1 + \cdots + f_{N-2} + \frac{f_{N-1}}{2}
and the constant contributes as $\int c = c(N-1)$, otherwise
.. math::
\int_0^L f(x) \d{x} \approx \sum_{n=0}^{N-1} f_n
and the constant constributes as $\int c = cN$.
"""
f_sum = sum(f)
if self.inclusive:
f_sum -= (f[0] + f[-1])/2
c = len(f) - 1
else:
c = len(f)
return f - f_sum/c
def R(self, f):
r"""Restriction operator:
.. math::
\mat{R} = \tfrac{1}{2}\mat{I}^{T} = \frac{1}{2}\begin{pmatrix}
1 & \tfrac{1}{2}\\
& \tfrac{1}{2} & 1 & \tfrac{1}{2}\\
&& \ddots\\
& & & & \tfrac{1}{2} & 1 & \tfrac{1}{2}\\
& & & & & & 1
\end{pmatrix}
"""
res = 1*f[::2]
res[:-1] += f[1::2]/2
res[1:] += f[1::2]/2
res /= 2
return res
def I(self, f):
r"""Interpolation operator:
.. math::
\mat{I} = \begin{pmatrix}
1\\
\tfrac{1}{2} & \tfrac{1}{2}\\
& 1 \\
& \tfrac{1}{2} & \tfrac{1}{2}\\
& & \ddots\\
& & & & \tfrac{1}{2} & \tfrac{1}{2}\\
& & & & & 1
\end{pmatrix}
"""
N = len(f)
N_new = 2*N - 1
res = np.empty(N_new, dtype=float)
res[::2] = f
res[1::2] = (f[1:] + f[:-1])/2
return res
def T(self, f, diag=False):
r"""Act with the operator $T(f)=f''$ at the highest level.
.. math::
\mat{T}_{\text{Neumann}} = \frac{1}{\delta^2}\begin{pmatrix}
-1 & 1 & \\
1 & -2 & 1\\
& \ddots & \ddots & \ddots\\
& & 1 & -2 & 1\\
& & & 1 & -1 \\
\end{pmatrix}.
.. note:: In order to preserve the Galerkin properties, we must
ensure that the step size $\delta_n = 2\delta_{n-1}$ even
though the lattice will not satisfy this except at the
highest level.
Parameters
----------
f : array
diag : bool
If `True`, then return the diagonals.
"""
m = int(self.get_N(len(f), inverse=True))
dx = self.dx_n/2**m
if diag:
res = -2*np.ones(len(f), dtype=f.dtype)
res[0] = -1
res[-1] = -1
else:
f_aug = np.empty(len(f)+2, dtype=f.dtype)
f_aug[1:-1] = f
f_aug[0] = f[0]
f_aug[-1] = f[-1]
res = -2*f + f_aug[2:] + f_aug[:-2]
res /= dx**2
return res
def D(self, f, x):
r"""Return $D(f)$, the diagonal portion of the problem (not including
the operator $-T = -\nabla^2$) acting on `f`.
This simply calls ``self.problem.D`` and then corrects the
endpoints if :attr:`inclusive`. It is intended that this only
be called on the top-level problem.
"""
d = self.problem.D(f, x)
if self.inclusive:
# Correct endpoints
d[[0, -1]] /= 2
return d
def smooth(self, f, b):
r"""Smoothing operator. The method depends on :attr:`smoothing_method`:
`0` : Weighted Jacobi:
.. math::
\vect{f} \rightarrow \vect{f} - \tfrac{2}{3}\mat{d}^{-1}
(\mat{A}\vect{f} - \vect{b})
`1` : Red-Black Gauss-Seidel:
If :attr:`inclusive` is `True`, then the abscissa
include the endpoints:
.. math::
x_{n} = \delta n, \qquad
\delta = \frac{L}{N-1}
::
x x x ... x x
|---|---|-- ... |---|
0 d 2d ... L
0 1 2 ...N-2 N-1
R---B---R-- ... B---R
If :attr:`inclusive` is `False`, then the abscissa are
interior:
.. math::
x_{n} = \frac{\delta}{2} + \delta n, \qquad
\delta = \frac{L}{N}
::
o x x x ... x x o
|--|----|----|-- ... |----|--|
0 d/2 3d/2 5d/2 ... L
0 1 2 ...N-2 N-1
---R----B----R-- ... B----R---
Implemented only for weight = 1
.. math::
\vect{z}^{m+1}_{red} = \mat{d}^{-1} (\vect{f}^{m} - A_{black} \vect{f}^{m}_{black}
)\\
\vect{z}^{m+1}_{black} = \mat{d}^{-1} (\vect{f}^{m} - A_{red} \vect{f}^{m}_{red}
)\\
\vect{z}^{m+1}_{black} = \mat{d}^{-1} (\vect{f}^{m} - A_{red} \vect{f}^{m}_{red}
)\\
\vect{z}^{m+1}_{red} = \mat{d}^{-1} (\vect{f}^{m} - A_{black} \vect{f}^{m}_{black}
)
`2` : Lexicographic Gauss-Seidel:
.. math::
\vect{z}^{m+1}_{red} = \mat{d}^{-1} (\vect{b} - \mat{L}\vect{f}^{m+1}
- \mat{U}\vect{f}^{m})\\
\vect{f}^{m+1} = \vect{f}^{m} + weight*(\vect{z}^{m+1}-\vect{f}^{m})
\vect{z}^{m+2}_{red} = \mat{d}^{-1} (\vect{b} - \mat{U}\vect{f}^{m+2}
- \mat{L}\vect{f}^{m+1})\\
\vect{f}^{m+2} = \vect{f}^{m+1} + weight*(\vect{z}^{m+2}-\vect{f}^{m+1})
Parameters
----------
f : array
Solution to be smoothed.
b : array
Right hand side.
"""
m = int(self.get_N(len(f), inverse=True))
res = -self.residue(f, b)
if 0 == self.smoothing_method:
# Weighted Jacobi
d = (-self.T(f, diag=True)
+ self.problem.D(0*f + 1, self.x(len(f))))
return f - 2/3*res/d
elif 1 == self.smoothing_method:
#weighted Red-Black Gauss-Seidel
#f[red] = d**(-1)(b[red]-A_black*f[black])
# = d**(-1)(b[red]+T_black*f[black])
# = d**(-1)(b[red]+f[black]/dx**2)
#and similarly for black
weight = 1
if weight != 1:
raise(NotImplementedError)
if ((f.shape[0]==1) |(np.mod(f.shape[0], 2)==0)):
raise(NotImplementedError)
dx = self.dx_n/2**m
d = (-self.T(f, diag=True)
+ self.problem.D(0*f + 1, self.x(len(f))))
f[0::2] = (b[0::2] + np.append(f[1::2], 0)/dx**2
+ np.append(0, f[1:-1:2])/dx**2)/d[0::2]
f[1::2] = (b[1::2] + f[0:-1:2]/dx**2 + f[2::2]/dx**2)/d[1::2]
f[1::2] = (b[1::2] + f[0:-1:2]/dx**2 + f[2::2]/dx**2)/d[1::2]
f[0::2] = (b[0::2] + np.append(f[1::2], 0)/dx**2
+ np.append(0, f[1:-1:2])/dx**2)/d[0::2]
#we operate twice to maintain symmetry of the operator.
return f
else:
#weighted Lexicographic Gauss-Seidel
weight = 1
d = (-self.T(f, diag=True)
+ self.problem.D(0*f + 1, self.x(len(f))))
dx = self.dx_n/2**m
f[0] = (1-weight)*f[0] + weight*(b[0]+f[1]/dx**2)/d[0]
for i in range(1, f.shape[0]-1):
f[i] = (1-weight)*f[i] +\
weight*(b[i]+f[i+1]/dx**2 + f[i-1]/dx**2)/d[i]
f[-1] = (1-weight)*f[-1] +\
weight*(b[-1] + f[-2]/dx**2)/d[-1]
return f
'''
.. plot::
:include-source:
import scipy.integrate
n0 = 2
n = 7
N = 2**n*(n0 - 1) + 1
Rs = {}
Is = {}
As = {}
one = ones(N, dtype=float)
T = np.diag(-2*one) + np.diag(one[:-1], 1) + np.diag(one[:-1], -1)
Neumann = False
if Neumann:
T[0, :2] = [-1, 1]
T[-1, -2:] = [1, -1]
L = 1
else:
T[0, -1] = 1
T[-1, 0] = 1
L = 2
dx = L/N
x = np.linspace(0, 1, N)*(L-dx) + dx/2
T /= dx**2
a = 5
f = np.exp(a*(cos(np.pi*x) - 1))
b = a*np.pi**2*(np.cos(np.pi*x) - a*np.sin(np.pi*x)**2)*f
f -= sp.integrate.quad(lambda x:np.exp(a*(cos(np.pi*x) - 1)),
0, L)[0]/L
def R(f):
res = 1*f[::2]
res[:-1] += f[1::2]/2
res[1:] += f[1::2]/2
res /= 2
return res
def I(f):
N = len(f)
N_new = 2*N - 1
res = np.empty(N_new, dtype=float)
res[::2] = f
res[1::2] = (f[1:] + f[:-1])/2
return res
def A(f):
"Apply the operator A"
n = 0
while len(f) < N:
f = I(f)
n += 1
f = np.dot(-T, f)
while n > 0:
f = R(f)
n -= 1
return f
As = []
Rs = []
Is = []
Ss = []
for m in xrange(n+1):
M = int(self.get_N(len(f)))
print m, n
zero = np.zeros(M, dtype=float)
Rs.append(np.zeros(((M+1)/2, M), dtype=float))
Is.append(np.zeros((2*M-1, M), dtype=float))
As.append(np.zeros((M, M), dtype=float))
for i_ in xrange(M):
zero[i_] = 1
Rs[-1][:, i_] = R(zero)
Is[-1][:, i_] = I(zero)
As[-1][:, i_] = A(zero)
zero[i_] = 0
Ss.append(np.eye(M)
- 2/3/diag(As[-1])[:, None]*As[-1])
bs = [b]
f_ = np.linalg.lstsq(As[-1], bs[0])[0]
f_ -= sum(f_)/len(f_)
fs = [f_]
for i_ in reversed(xrange(n)):
bs.insert(0, R(bs[0]))
m = int(self.get_N(len(f), inverse=True))
fs.insert(0, np.linalg.lstsq(As[m], bs[0])[0])
fs[0] -= sum(fs[0])/len(fs[0])
for i_ in xrange(n+1):
ans = np.linalg.lstsq(As[i_], bs[i_])[0]
while len(ans) < N:
ans = I(ans)
ans -= sum(ans)/len(ans)
plot(x, ans - fs[-1])
def smooth(f, b):
m = int(self.get_N(len(f), inverse=True))
res = A(f) - b
return f - 2/3*res/np.diag(As[m])
def v_cycle(f, b):
m = int(self.get_N(len(f), inverse=True))
if m == 0:
A_ = As[m]
f = np.linalg.lstsq(A_, b)[0]
else:
f = smooth(f, b)
f = smooth(f, b)
err = A(f) - b
corr = I(v_cycle(0*R(f), R(err)))
f -= corr
f = smooth(f, b)
f = smooth(f, b)
return f - sum(f)/len(f)
def fmg():
f = v_cycle(0*bs[0], bs[0])
for m in xrange(1, n+1):
f = v_cycle(I(f), bs[m])
return f
'''
def comparing_periodic():
'''comparing the various implementation of periodic'''
mg = MultigridPeriodic(inclusive=False, smoothing_method=0, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(121)
pl.plt.semilogy(res, 'r:', label='residue J')
pl.plt.semilogy(err_, 'r-', label='err J')
pl.plt.semilogy(err, 'r--', label='exact err J')
pl.plt.legend()
[docs]def comparing_periodic():
'''comparing the various smoothings of periodic'''
mg = MultigridPeriodic(inclusive=False, smoothing_method=0, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.semilogy(res, 'r:', label='residue J')
pl.plt.semilogy(err_, 'r-', label='err J')
pl.plt.semilogy(err, 'r--', label='exact err J')
pl.plt.legend()
mg = MultigridPeriodic(inclusive=False, smoothing_method=1, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.semilogy(res, 'g:', label='residue RBGS')
pl.plt.semilogy(err_, 'g-', label='err RBGS')
pl.plt.semilogy(err, 'g--', label='exact err RBGS')
pl.plt.legend()
mg = MultigridPeriodic(inclusive=False, smoothing_method=2, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.semilogy(res, 'b:', label='residue LGS')
pl.plt.semilogy(err_, 'b-', label='err LGS')
pl.plt.semilogy(err, 'b--', label='exact err LGS')
pl.plt.legend()
pl.plt.title('Not Inclusive')
def comparing_neumann():
'''comparing the various implementation of neumann'''
mg = MultigridNeumann(inclusive=False, smoothing_method=0, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(121)
pl.plt.semilogy(res, 'r:', label='residue J')
pl.plt.semilogy(err_, 'r-', label='err J')
pl.plt.semilogy(err, 'r--', label='exact err J')
pl.plt.legend()
mg = MultigridNeumann(inclusive=False, smoothing_method=1, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(121)
pl.plt.semilogy(res, 'g:', label='residue RBGS')
pl.plt.semilogy(err_, 'g-', label='err RBGS')
pl.plt.semilogy(err, 'g--', label='exact err RBGS')
pl.plt.legend()
mg = MultigridNeumann(inclusive=False, smoothing_method=2, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(121)
pl.plt.semilogy(res, 'b:', label='residue LGS')
pl.plt.semilogy(err_, 'b-', label='err LGS')
pl.plt.semilogy(err, 'b--', label='exact err LGS')
pl.plt.legend()
pl.plt.title('Not Inclusive')
mg = MultigridNeumann(inclusive=True, smoothing_method=0, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
b[[0, -1]] /= 2
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(122)
pl.plt.semilogy(res, 'r:', label='residue J')
pl.plt.semilogy(err_, 'r-', label='err J')
pl.plt.semilogy(err, 'r--', label='exact err J')
pl.plt.legend()
mg = MultigridNeumann(inclusive=True, smoothing_method=1, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
b[[0, -1]] /= 2
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(122)
pl.plt.semilogy(res, 'g:', label='residue RBGS')
pl.plt.semilogy(err_, 'g-', label='err RBGS')
pl.plt.semilogy(err, 'g--', label='exact err RBGS')
pl.plt.legend()
pl.plt.title('Inclusive')
mg = MultigridNeumann(inclusive=True, smoothing_method=2, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
b[[0, -1]] /= 2
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*x
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(10):
g = mg.f_cycle(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
pl.plt.subplot(122)
pl.plt.semilogy(res, 'b:', label='residue LGS')
pl.plt.semilogy(err_, 'b-', label='err LGS')
pl.plt.semilogy(err, 'b--', label='exact err LGS')
pl.plt.legend()
pl.plt.title('Inclusive')
def test_smoother_convergence(smoothing_method = 0, inclusive=False, periodic=True):
''' Testing the convergence of various smoothers for inclusive and non-inclusive boundary
conditions '''
if (periodic):
mg = MultigridPeriodic(inclusive=False, smoothing_method=smoothing_method, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*f
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(1000):
g = mg.smooth(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
else:
mg = MultigridNeumann(inclusive=inclusive, smoothing_method=smoothing_method, neutralize=True)
x = mg.x(); dx = x[1]-x[0]
f = mg.problem.f(x)
b = mg.problem.b(x)
A = mg.to_mat(mg.A, len(f))
if inclusive:
b[[0, -1]] /= 2
f_ = mg.make_neutral(np.linalg.lstsq(A, b)[0])
g = 0*f
res = [abs(mg.residue(g, b)).max()]
err = [abs(g - f).max()]
err_ = [abs(g - f_).max()]
for n in xrange(1000):
g = mg.smooth(g, b)
res.append(abs(mg.residue(g, b)).max())
err.append(abs(g - f).max())
err_.append(abs(g - f_).max())
#pl.clf()
pl.subplot(211)
pl.plot(res, ':', label='residue')
pl.subplot(212)
pl.plot(err_, '-', label='err')
pl.plot(err, '--', label='exact err')
pl.legend()
pl.title('Iterating using smoother', )
def test_I_R(m=3, periodic=False, inclusive=False):
""" testing the interpolation and the restriction operator for various boundary conditions """
n = m+1
if periodic:
mg = MultigridPeriodic(inclusive=False, n=n)
else:
mg = MultigridNeumann(inclusive=False, n=n)
mats = mg.mats()
Rmn = mats.R[n]
An, Am = mats.A[n], mats.A[m]
Inm = mats.I[m]
print(Am - np.dot(np.dot(Rmn, An), Inm))
class Neumann(object):
"""Demo of 1D Neumann problem."""