Source code for dec.spectral

'''
Spectral DEC
=============
'''
from dec.helper import *
from numpy import *
from dec.forms import *
from numpy.testing import assert_almost_equal
import operator
from functools import reduce
try:
    from scipy.linalg.special_matrices import toeplitz
except ImportError:
    from scipy.linalg.basic import toeplitz

seterr(divide='ignore', invalid='ignore')

[docs]def alpha0(N, x): r''' .. math:: \alpha_{N}(x)=\frac{1}{N} \begin{cases} \cot\frac{x}{2}\,\sin\frac{Nx}{2} & \text{if N even,}\\ \csc\frac{x}{2}\,\sin\frac{Nx}{2} & \text{if N odd.} \end{cases} >>> def a0(N, i): return round(alpha0(N, i*2*pi/N), 15) >>> (a0(5, 0), a0(5, 1), a0(5, 2)) == (1.0, 0.0, 0.0) True >>> (a0(6, 0), a0(6, 1), a0(6, 2)) == (1.0, 0.0, 0.0) True ''' if N % 2 == 0: y = (sin(N*x/2) / tan(x/2)) / N else: y = (sin(N*x/2) / sin(x/2)) / N if hasattr(y, '__setitem__'): y[x==0] = 1 else: if x==0: y = 1 return y
[docs]def beta0(N, x): r''' .. math:: \beta_{N}(x)= \begin{cases} \frac{1}{2\pi} -\frac{1}{4}\cos\frac{Nx}{2} + \frac{1}{N}\sum\limits_{n=1}^{N/2} \frac{n\cos nx}{\sin\frac{n\pi}{N}} & \text{if N even,}\\ \frac{1}{2\pi} + \frac{1}{N}\sum\limits_{n=1}^{(N-1)/2} \frac{n\cos nx}{\sin\frac{n\pi}{N}} & \text{if N odd.} \end{cases} ''' if N % 2 == 0: y = 1/(2*pi) - cos(N*x/2)/4 for n in range(1, N//2 + 1): y += n*cos(n*x) / sin(n*pi/N) / N else: y = 1/(2*pi) + 0*x for n in range(1, (N-1)//2 + 1): y += n*cos(n*x) / sin(n*pi/N) / N return y
[docs]def alpha(N, n, x): r''' .. math:: \alpha_{N, n} (x) = \alpha_N(x-h n) ''' return alpha0(N, x-2*pi/N*n)
[docs]def beta(N, n, x): r''' .. math:: \beta_{N, n} (x) = \beta_N(x-h n) ''' return beta0(N, x-2*pi/N*n)
[docs]def gamma(N, k): r''' .. math:: \gamma_{N, n} = \int_0^{\frac{h}{2}=\frac{\pi}{N}} \beta_{N,n}(x) \, dx = \begin{cases} \frac{1}{2N} - \frac{(-1)^k}{2N} + \frac{1}{N}\sum\limits_{n=1}^{N/2} \frac{\sin(2 k n \pi/N) - \sin((2 k-1) n \pi/N)} {\sin\frac{n\pi}{N}} & \text{if N even,}\\ \frac{1}{2N} + \frac{1}{N}\sum\limits_{n=1}^{(N-1)/2} \frac{\sin(2 k n \pi/N) - \sin((2 k-1) n \pi/N)}{\sin\frac{n\pi}{N}} & \text{if N odd.} \end{cases} ''' if N % 2 == 0: y = 1/(2*N) - (-1)**k/(2*N) for n in range(1, N//2 + 1): y += (sin(2*n*k*pi/N) - sin((2*k-1)*n*pi/N)) / sin(n*pi/N) / N else: y = 1/(2*N) for n in range(1, (N-1)//2 + 1): y += (sin(2*n*k*pi/N) - sin((2*k-1)*n*pi/N)) / sin(n*pi/N) / N return y ######################################## # Periodic Basis Functions ########################################
[docs]def phi0(N, n, x): r''' Basis functions for primal 0-forms. .. math:: \phi_{N,n}^{0}(x)=\alpha_{N,n}(x) ''' return alpha(N, n, x)
[docs]def phi1(N, n, x): r''' Basis functions for primal 1-forms. .. math:: \phi_{N,n}^{1}(x)=\beta_{N,n+\frac{1}{2}}(x) ''' return beta(N, n + 0.5, x)
[docs]def phid0(N, n, x): r''' Basis functions for dual 0-forms. .. math:: \widetilde{\phi}_{N,n}^{0}(x)=\alpha_{N,n+\frac{1}{2}}(x) ''' return alpha(N, n + 0.5, x)
[docs]def phid1(N, n, x): r''' Basis functions for dual 1-forms. .. math:: \widetilde{\phi}_{N,n}^{1}(x)=\alpha_{N,n}(x) ''' return beta(N, n, x) ######################################## # Mapping between semi-circle and line # ########################################
[docs]def varphi(x): r''' .. math:: \varphi:&& [-1,1]\to[0,\pi]\\ && x\mapsto\arccos(-x) >>> varphi(-1) == 0 True >>> varphi(+1) == pi True ''' return arccos(-x)
[docs]def varphi_inv(x): r''' .. math:: \varphi^{-1}:&& [0,\pi]\to[-1,1]\\ && \theta\mapsto-\cos(\theta) >>> varphi_inv(0) == -1 True >>> varphi_inv(pi) == +1 True ''' return -cos(x) ############################################################ # Regular Basis Functions ############################################################
[docs]def delta(N, x): r''' Basis function for dual half-edge at the boundary. .. math:: \delta_{N}(x)=\left((N-1)^{2}\alpha_{2N-2,0}(x)+\frac{1}{2}\cos\left((N-1)x\right)\right)\sin(x) ''' return ((N-1)**2 * alpha(2*N-2, 0, x) + 0.5*cos((N-1)*x))*sin(x)
[docs]def rho(N, n, x): r''' .. math:: \rho_{N, n}(x) = 2 \gamma_{2N-2, n}\delta_N(x) + 2 \gamma_{2N-2, N-n-1} \delta_N(\pi-x) ''' return 2*gamma(2*N-2, n)*delta(N, x) + 2*gamma(2*N-2, N-n-1)*delta(N, pi-x)
[docs]def kappa0(N, n, x): r''' Basis functions for primal 0-forms. .. math:: \kappa_{N,n}^{0}(x) = \begin{cases} \alpha_{2N-2,n}(x), & n\in\{0,N-1\}\\ \alpha_{2N-2,n}(x)+\alpha_{2N-2,2N-2-n}(x), & n\in\{1,\dots,N-2\} \end{cases} ''' if n in (0, N-1): return alpha(2*N-2, n, x) else: return (alpha(2*N-2, n, x) + alpha(2*N-2, 2*N-2-n, x))
[docs]def kappa1(N, n, x): r''' .. math:: \kappa_{N,n}^{1}(x) = \left( \beta_{2N-2,n+\frac{1}{2}}(x)- \beta_{2N-2,2N-3-n+\frac{1}{2}}(x)\right) \mathbf{d}x,n\in\{0,\dots,N-2\} ''' return (beta(2*N-2, n+0.5, x) - beta(2*N-2, 2*N-3-n+0.5, x))
[docs]def kappad0(N, n, x): r''' .. math:: \widetilde{\kappa}_{N,n}^{0}(x)= \alpha_{2N-2,\, n+\frac{1}{2}}(x)+ \alpha_{2N-2,\,2N-3-n+\frac{1}{2}}(x),\quad n\in\{0,\dots,N-2\} ''' return (alpha(2*N-2, n+0.5, x) + alpha(2*N-2, 2*N-3-n+0.5, x))
[docs]def kappad1(N, n, x): r''' .. math:: \widetilde{\kappa}_{N,n}^{1}(x)= \begin{cases} \delta(x)\mathbf{d}x & \qquad n=0\\ \left(\beta_{2N-2,n}(x)-\beta_{2N-2,2N-2-n}(x)-\rho_{N,n}(x)\right)\mathbf{d}x & \qquad n\in\{1,\dots,N-2\}\\ \delta(\pi-x)\mathbf{d}x & \qquad n=N-1 \end{cases} ''' if n == 0: return delta(N, x) if n == N-1: return delta(N, pi - x) y = (beta(2*N-2, n, x) - beta(2*N-2, 2*N-2-n, x) - rho(N, n, x) ) return y ############################################################ # Chebyshev Basis Functions # They seem to be all polynomials: # >>> x = linspace(-1, +1, 10000) # >>> round_(polyfit(x, psid1(5, 5, x), 10), 3)[::-1] # array([ 0.125, 0.657, -2.177, -9.941, -5.121, 9.941, 7.877, -0. , # -0. , 0. , 0. ]) ############################################################ #TODO: This is a dirty hack. Compute the limits explicilty?
def __fix_singularity_at_boundary(x): # Avoid division by zero # Make sure function has no side effects if type(x) is ndarray: x += (x==-1)*1e-16 - (x==+1)*1e-16 else: if x==-1: x += 1e-16 if x==+1: x -= 1e-16 return x
[docs]def psi0(N, n, x): r''' Basis functions for primal 0-forms. .. math:: \psi_{N,n}^{0}(x)=\kappa_{N,n}^{0}(\arccos(-x)) ''' return kappa0(N, n, arccos(-x))
[docs]def psi1(N, n, x): r''' Basis functions for primal 1-forms. .. math:: \psi_{N,n}^{1}(x)\mathbf{d}x= \kappa_{N,n}^{1}(\arccos(-x))\frac{\mathbf{d}x}{\sqrt{1-x^{2}}} ''' #x = __fix_singularity_at_boundary(x) return kappa1(N, n, arccos(-x))/sqrt(1 - x**2)
[docs]def psid0(N, n, x): r''' Basis functions for dual 0-forms. .. math:: \tilde{\psi}_{N,n}^{0}(x)=\tilde{\kappa}_{N,n}^{0}(\arccos(-x)) ''' return kappad0(N, n, arccos(-x))
[docs]def psid1(N, n, x): r''' Basis functions for dual 1-forms. .. math:: \tilde{\psi}_{N,n}^{1}(x)\mathbf{d}x=\tilde{\kappa}_{N,n}^{1}(\arccos(-x))\frac{\mathbf{d}x}{\sqrt{1-x^{2}}} ''' x = __fix_singularity_at_boundary(x) return kappad1(N, n, arccos(-x))/sqrt(1 - x**2) ########################### # Lagrange polynomials ###########################
[docs]def lagrange_polynomials(x): r''' Lagrange Polynomials for the set of points defined by :math:`x_m`. The Lagrange Polynomials are such that they are 1 at the point, and 0 everywhere else. .. math:: \psi_{n}^{0}(x)=\prod_{m=0,m\neq n}^{N-1}\frac{x-x_{m}}{x_{n}-x_{m}} >>> L = lagrange_polynomials([0, 1, 2]) >>> [l(0) for l in L] [1.0, 0.0, -0.0] >>> [l(1) for l in L] [-0.0, 1.0, 0.0] >>> [l(2) for l in L] [0.0, -0.0, 1.0] ''' N = len(x) L = [None]*N for i in range(N): def Li(y, i=i): gen = ((y-x[j])/(x[i]-x[j]) for j in set(range(N)).difference([i])) return reduce(operator.mul, gen) L[i] = (Li) return L
def freq(N): ''' >>> freq(5) array([-2., -1., 0., 1., 2.]) >>> freq(6) array([-3., -2., -1., 0., 1., 2.]) ''' return fft.fftshift(fft.fftfreq(N)*N) F = lambda x: fft.fftshift(fft.fft(x)) Finv = lambda x: fft.ifft(fft.ifftshift(x)) def method_in_fourier_space(g): def f(x, *args, **keywords): return Finv(g(F(x), *args, **keywords)) f.__doc__ = g.__doc__ # make sure doctests are run return f
[docs]def H(a): r''' .. math:: \mathbf{H}^0 = \widetilde{\mathbf{H}}^0 = \mathcal{F}^{-1} \mathbf{I}^{-\frac{h}{2}, \frac{h}{2}} \mathcal{F} ''' N = len(a) h = 2*pi/N return Finv( F(a) * I_diag(N, -h/2, h/2) )
[docs]def Hinv(a): r''' .. math:: \mathbf{H}^1 = \widetilde{\mathbf{H}}^1 = \mathcal{F}^{-1}{\mathbf{I}^{-\frac{h}{2},\frac{h}{2}}}^{-1}\mathcal{F} ''' N = len(a) h = 2*pi/N return Finv( F(a) / I_diag(N, -h/2, h/2) )
def I(a, x0, x1): N = len(a) return Finv ( I_diag(N, x0, x1) * F(a) ) def Iinv(a, x0, x1): N = len(a) return Finv ( F(a) / I_diag(N, x0, x1) ) def S(a): N = len(a) h = 2*pi/N return Finv( F(a) * S_diag(N, h/2) ) def Sinv(a): N = len(a) h = 2*pi/N return Finv( F(a) * S_diag(N, -h/2) )
[docs]def I_diag(N, a, b): r''' A diagonal matrix that corresponds to integration in Fourier space. Corresponds to :math:`f(x) \mapsto \int_{x+a}^{x+b} f(\xi) d\xi` .. math:: \mathbf{I}_{\phantom{a,b}nn}^{a,b} =\frac{e^{inb}-e^{ina}}{in} .. math:: \mathbf{I}_{\phantom{a,b}00}^{a,b}=b-a ''' n = freq(N) y = (exp(1j*n*b) - exp(1j*n*a))/(1j*n) y[n==0] = b - a return y
[docs]def S_diag(N, a): r''' A diagonal matrix that corresponds to shifting in Fourier Space Corresponds to :math:`f(x) \mapsto f(x-h)` .. math:: \mathbf{S}_{\phantom{a}nn}^{a}=e^{ina} ''' n = freq(N) return exp(1j*n*a)
def fourier_I(x, a, b): r''' Corresponds to :math:`f(x) \mapsto \int_{x+a}^{x+b} f(\xi) d\xi` ''' N = x.shape[0] n = freq(N) y = (exp(1j*n*b) - exp(1j*n*a))/(1j*n) y[n==0] = b - a return x*y def fourier_I_inv(x, a, b): N = x.shape[0] n = freq(N) y = (exp(1j*n*b) - exp(1j*n*a))/(1j*n) y[n==0] = b - a return x/y def fourier_T(x, h): r''' Corresponds to :math:`f(x) \mapsto f(x+a)` ''' N = x.shape[0] n = freq(N) return x*exp(1j*n*h) def fourier_J(x, a, b, c): r''' Corresponds to :math:`f(x) \mapsto \int_{x+a}^{x+b} f(\xi) \sin(\xi+c) d\xi` ''' c = fourier_I(x, a, b) b = (roll(c[2:]*x,+1)*exp(1j*c) - roll(c[:-2]*x,-1)*exp(-1j*c))/2j return b def refine(x): ''' Resample x at a twice refined grid. >>> N = 4 >>> x = linspace(0, 2*pi, N+1)[:-1] >>> y = linspace(0, 2*pi, 2*N+1)[:-1] >>> approx(refine(cos(x)), cos(y)) True ''' x = interweave(x, S(x)) return x def subdivide(x): ''' Subdivide x like below. Assume points are equidistant. * * * * * * * * * * * * * * * * * ''' if len(x) < 2: return x assert is_equidistant(x) y = interweave(x[:-1], x[:-1] + 0.5*diff(x)) return concatenate((y, [x[-1]])) def integrate_spectral_coarse(x, f): ''' >>> integrate_spectral_coarse(linspace(-pi, pi, 3), sin) array([-2., 2.]) ''' assert is_equidistant(x) assert approx(2*pi, x[-1] - x[0]), x f0 = f(x[:-1] + 0.5*diff(x)) f1 = real(H(f0)) return f1 def integrate_spectral(x, f): ''' >>> integrate_spectral(linspace(-pi, pi, 3), sin) array([-2., 2.]) ''' assert is_equidistant(x) r = subdivide f1 = integrate_spectral_coarse(r(r(r(x))), f) return (f1[0::8] + f1[1::8] + f1[2::8] + f1[3::8] + f1[4::8] + f1[5::8] + f1[6::8] + f1[7::8]) def integrate_chebyshev(xi, f): ''' >>> integrate_chebyshev(array([cos(pi), cos(.5*pi), cos(0)]), lambda x: x) array([-0.5, 0.5]) ''' assert approx(2, xi[-1] - xi[0]), xi F = lambda theta: f(-cos(theta))*sin(theta) x = varphi(xi) assert(is_equidistant(x)) # complete the circle from the other side x = concatenate((x, x[1:]+pi)) return integrate_spectral(x, F)[:len(xi)-1] def integrate_chebyshev_dual(xi, f): ''' Integrate points that may include the two half-edges at the boundary. #>>> integrate_chebyshev_dual(array([cos(pi), cos(0.75*pi), cos(0.25*pi), cos(0)]), lambda x: x) #array([-0.25, 0. , 0.25]) ''' x = varphi(xi) z = varphi_inv(concatenate(([x[0]], subdivide(x[1:-1]), [x[-1]]))) i, j = concatenate(( [0], integrate_chebyshev(z, f), [0] )).reshape(-1, 2).T return i+j def split_args(I): ''' Convert integration function from I(x, f) to I(x0, x1, f) form ''' def J(x0, x1, f, I=I): assert_almost_equal(x0[1:], x1[:-1]) return I(concatenate((x0, [x1[-1]])), f) return J ############################## # Discrete exterior calculus ############################## def hodge_star_matrix(P, B): ''' Compute Hodge-Star matrix from basis functions. ''' P0, P1, P0d, P1d = P B0, B1, B0d, B1d = B H0 = vstack(P1d(b) for b in B0).T H1 = vstack(P0d(b) for b in B1).T H0d = vstack(P1(b) for b in B0d).T H1d = vstack(P0(b) for b in B1d).T return H0, H1, H0d, H1d def reconstruction(basis_fn): ''' Give the reconstruction functions for the set of basis functions basis_fn. ''' def R(a, B): def r(*x): return sum(a[i]*f(*x) for i, f in enumerate(B)) return r return [(lambda a, B=B: R(a, B)) for B in basis_fn] Grid_1D_Interface = ''' dimension xmin xmax delta delta_dual verts verts_dual edges edges_dual basis_fn projection derivative hodge_star '''.split() def check(g, interface): ''' Check whether an object satisfies an interface. >>> check(Grid_1D_Periodic(4), Grid_1D_Interface) True >>> check(Grid_1D_Chebyshev(4), Grid_1D_Interface) True >>> check(Grid_1D_Regular(4), Grid_1D_Interface) True ''' rslt = True for i in interface: rslt = (rslt and hasattr(g, i)) return rslt
[docs]class Grid_1D_Periodic(object): def __init__(self, n, xmin=0, xmax=2*pi): assert xmax > xmin dimension = 1 pnts = linspace(xmin, xmax, num=(n+1)) lenx = abs(xmax - xmin) delta = diff(pnts) verts = pnts[:-1] edges = (pnts[:-1], pnts[1:]) verts_dual = verts + 0.5*delta edges_dual = (roll(verts_dual,shift=1), verts_dual) edges_dual[0][0] -= lenx delta_dual = delta V = verts S0 = arange(len(V)) S1 = (S0[:-1], S0[1:]) Vd = verts_dual S0d = arange(len(Vd)) S1d = (S0d[:-1], S0d[1:]) self.dimension = dimension self.n = n self.xmin = xmin self.xmax = xmax self.delta = delta self.delta_dual = delta_dual self.pnts = pnts self.verts = verts self.verts_dual = verts_dual self.edges = edges self.edges_dual = edges_dual def projection(self): P0 = lambda f: f(self.verts) P1 = lambda f: split_args(integrate_spectral)( self.edges[0], self.edges[1], f) P0d = lambda f: f(self.verts_dual) P1d = lambda f: split_args(integrate_spectral)( self.edges_dual[0], self.edges_dual[1], f) return P0, P1, P0d, P1d def basis_fn(self): n = self.n B0 = [lambda x, i=i: phi0(n, i, x) for i in range(n)] B1 = [lambda x, i=i: phi1(n, i, x) for i in range(n)] B0d = [lambda x, i=i: phid0(n, i, x) for i in range(n)] B1d = [lambda x, i=i: phid1(n, i, x) for i in range(n)] return B0, B1, B0d, B1d def reconstruction(self): R0, R1, R0d, R1d = reconstruction(self.basis_fn()) return R0, R1, R0d, R1d def derivative(self): D0 = lambda f: roll(f, shift=-1) - f D0d = lambda f: roll(D0(f), shift=+1) return D0, D0d def hodge_star(self): H0 = lambda x: real(H(x)) H1 = lambda x: real(Hinv(x)) H0d = H0 H1d = H1 return H0, H1, H0d, H1d def derivative_matrix(self): n = self.n rng = arange(n) ons = ones(n) d = row_stack(( column_stack(( rng, roll(rng, shift= -1), +ons)), column_stack(( rng, rng, -ons)) )) D = sparse_matrix(d, n, n) return D, -D.T def differentiation_toeplitz(self): raise NotImplemented #TODO: fix this implementation n = self.n h = 2*pi/n assert n % 2 == 0 column = concatenate(( [ 0 ], (-1)**arange(1,n) / tan(arange(1,n)*h/2) )) row = concatenate(( column[:1], column[1:][::-1] )) D = toeplitz(column, row) return D def hodge_star_toeplitz(self): ''' The Hodge-Star using a Toeplitz matrix. ''' P0, P1, P0d, P1d = self.projection() column = P1d(lambda x: alpha0(self.n, x)) row = concatenate((column[:1], column[1:][::-1])) return toeplitz(column, row) def wedge(self): ''' Return \alpha ^ \beta. Keep only for primal forms for now. ''' def w00(alpha, beta): return alpha * beta def _w01(alpha, beta): return S(H( alpha * Hinv(Sinv(beta)) )) def w01(alpha, beta): # a = interweave(alpha, T(alpha, [S])) # b = interweave(T(beta, [Hinv, Sinv]), T(beta, [Hinv])) a = refine(alpha) b = refine(Hinv(Sinv(beta))) c = S(H(a * b)) return c[0::2] + c[1::2] return w00, w01, _w01 def contraction(self, V): ''' Return i_V. Keep only for primal forms for now. ''' def C1(alpha): return Hinv(Sinv(V)) * Hinv(Sinv(alpha)) return C1
def _slow_integration(a, b, f): from scipy.integrate import quad return array([quad(f, _a, _b)[0] for _a, _b in zip(a, b)])
[docs]def A_diag(N): r''' .. math:: \mathbf{A}=\text{diag}\left(\begin{array}{ccccccc}\frac{1}{2} & 1 & 1 & \dots & 1 & 1 & \frac{1}{2}\end{array}\right) ''' d = concatenate(([0.5], ones(N-2), [0.5])) return d
[docs]def H0_regular(f): r''' .. math:: \mathbf{H}^{0}= \mathbf{A} \mathbf{M}_{0}^{\dagger} \mathbf{I}^{-\frac{h}{2},\frac{h}{2}} \mathbf{M}_{0}^{+} ''' f = mirror0(f, +1) N = f.shape[0]; h = 2*pi/N f = I_space(-h/2, h/2)(f) f = unmirror0(f) f = f*A_diag(f.shape[0]) return real(f)
[docs]def H1_regular(f): r''' .. math:: \mathbf{H}^{1}= \mathbf{M}_{1}^{\dagger} {\mathbf{I}^{-\frac{h}{2},\frac{h}{2}}}^{-1} \mathbf{M}_{1}^{-} ''' f = mirror1(f, -1) N = f.shape[0]; h = 2*pi/N f = I_space_inv(-h/2, h/2)(f) f = unmirror1(f) return f
[docs]def H0d_regular(f): r''' .. math:: \widetilde{\mathbf{H}}^{0}= \mathbf{M}_{1}^{\dagger} \mathbf{I}^{-\frac{h}{2},\frac{h}{2}} \mathbf{M}_{1}^{-} ''' f = mirror1(f, -1) N = f.shape[0]; h = 2*pi/N f = I_space(-h/2, h/2)(f) f = unmirror1(f) return f
[docs]def H1d_regular(f): r''' .. math:: \widetilde{\mathbf{H}}^{1}= \mathbf{M}_{1}^{\dagger} {\mathbf{I}^{-\frac{h}{2},\frac{h}{2}}}^{-1} \mathbf{M}_{1}^{+} \mathbf{A}^{-1} ''' f = f/A_diag(f.shape[0]) f = mirror0(f, +1) N = f.shape[0]; h = 2*pi/N f = I_space_inv(-h/2, +h/2)(f) f = unmirror0(f) return f
[docs]class Grid_1D_Regular: def __init__(self, n, xmin=0, xmax=pi): assert xmax > xmin dimension = 1 N = n pnts = linspace(xmin, xmax, num=n) lenx = abs(xmax - xmin) delta = diff(pnts) verts = pnts verts_dual = verts[:-1] + 0.5*delta edges = (pnts[:-1], pnts[1:]) tmp = concatenate(([xmin], verts_dual, [xmax])) delta_dual = diff(tmp) edges_dual = (tmp[:-1], tmp[1:]) self.dimension = dimension self.n = n self.xmin = xmin self.xmax = xmax self.delta = delta self.delta_dual = delta_dual self.pnts = pnts self.verts = verts self.verts_dual = verts_dual self.edges = edges self.edges_dual = edges_dual def projection(self): P0 = lambda f: f(self.verts) P1 = lambda f: _slow_integration(self.edges[0], self.edges[1], f) P0d = lambda f: f(self.verts_dual) P1d = lambda f: _slow_integration(self.edges_dual[0], self.edges_dual[1], f) return P0, P1, P0d, P1d def basis_fn(self): n = self.n B0 = [lambda x, i=i: kappa0(n, i, x) for i in range(n)] B1 = [lambda x, i=i: kappa1(n, i, x) for i in range(n-1)] B0d = [lambda x, i=i: kappad0(n, i, x) for i in range(n-1)] B1d = [lambda x, i=i: kappad1(n, i, x) for i in range(n)] return B0, B1, B0d, B1d def reconstruction(self): R0, R1, R0d, R1d = reconstruction(self.basis_fn()) return R0, R1, R0d, R1d def derivative(self): D0 = lambda f: diff(f) D0d = lambda f: diff(concatenate(([0], f, [0]))) return D0, D0d def hodge_star(self): H0 = H0_regular H1 = H1_regular H0d = H0d_regular H1d = H1d_regular return H0, H1, H0d, H1d
[docs]def extend(f, n): r''' .. math:: \mathbf{E}^{n}:\quad \begin{bmatrix}x_{0}\\ \vdots\\ x_{N-1} \end{bmatrix} \mapsto \frac{N+2n}{N} \begin{bmatrix}\left.\begin{array}{c} 0\\ \vdots\\ 0 \end{array}\right\} n\\ \begin{array}{c} x_{0}\\ \vdots\\ x_{N-1} \end{array}\\ \left.\begin{array}{c} 0\\ \vdots\\ 0 \end{array}\right\} n \end{bmatrix} ''' N = f.shape[0] return (N+2.0*n)/N*_extend(f, n)
[docs]def unextend(f, n): r''' .. math:: \mathbf{E}^{-n}:\quad\begin{bmatrix}\left.\begin{array}{c} x_{-n}\\ \vdots\\ x_{-1} \end{array}\right\} n\\ \begin{array}{c} x_{0}\\ \vdots\\ x_{N-1} \end{array}\\ \left.\begin{array}{c} x_{N}\\ \vdots\\ x_{N+n-1} \end{array}\right\} n \end{bmatrix}\mapsto\frac{N-2n}{N}\begin{bmatrix}\left.\begin{array}{c} x_{0}+x_{N}\\ \vdots\\ x_{n-1}+x_{N+n-1} \end{array}\right\} n\\ \begin{array}{c} x_{n}\\ \vdots\\ x_{N-n-1} \end{array}\\ \left.\begin{array}{c} x_{N-n}+x_{-n}\\ \vdots\\ x_{N-1}+x_{-1} \end{array}\right\} n \end{bmatrix} ''' N = f.shape[0] return (N-2.0*n)/N*_unextend(f, n)
def _extend(f, n): ''' >>> _extend([1, 2, 3, 4], 2) array([ 0., 0., 1., 2., 3., 4., 0., 0.]) ''' return concatenate(( zeros(n), f, zeros(n) )) def _unextend(f, n): ''' >>> _unextend(array([0,0,1,2,3,4,0,0]), 2) array([1, 2, 3, 4]) >>> _unextend(array([1,2,0,0,3,4]), 2) array([4, 6]) >>> _unextend(array([0,1,2,3,4,5,6,7]), 2) array([ 8, 10, 4, 6]) ''' assert 2*n <= f.shape[0] x = f[n:-n].copy() x[:n] += f[-n:] x[-n:] += f[:n] return x
[docs]def mirror0(f, sign=+1): r''' .. math:: \mathbf{M}_{0}^{\pm}:\quad\begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-2}\\ x_{N-1} \end{bmatrix}\mapsto \begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-2}\\ x_{N-1}\\ \pm x_{N-2}\\ \vdots\\ \pm x_{1} \end{bmatrix} >>> mirror0(array([1, 2, 3])) array([1, 2, 3, 2]) >>> mirror0(array([1, 2, 3]), -1) array([ 1, 2, 3, -2]) >>> to_matrix(mirror0, 3) array([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.], [ 0., 1., 0.]]) >>> to_matrix(lambda x: mirror0(x, -1), 3) array([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.], [-0., -1., -0.]]) ''' return concatenate((f, sign*f[::-1][1:-1]))
[docs]def mirror1(f, sign=+1): r''' .. math:: \mathbf{M}_{1}^{\pm}:\quad \begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-2}\\ x_{N-1} \end{bmatrix}\mapsto \begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-2}\\ x_{N-1}\\ \pm x_{N-1}\\ \pm x_{N-2}\\ \vdots\\ \pm x_{1}\\ \pm x_{0} \end{bmatrix} >>> mirror1(array([1, 2, 3])) array([1, 2, 3, 3, 2, 1]) >>> mirror1(array([1, 2, 3]), -1) array([ 1, 2, 3, -3, -2, -1]) >>> to_matrix(mirror1, 2) array([[ 1., 0.], [ 0., 1.], [ 0., 1.], [ 1., 0.]]) >>> to_matrix(lambda x: mirror1(x, -1), 2) array([[ 1., 0.], [ 0., 1.], [-0., -1.], [-1., -0.]]) ''' return concatenate((f, sign*f[::-1]))
[docs]def unmirror0(f): r''' .. math:: \mathbf{M}_{0}^{\dagger}:\quad\begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-1}\\ x_{N}\\ \vdots\\ x_{2N-1} \end{bmatrix}\mapsto\begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-1}\\ x_{N}\\ x_{N+1} \end{bmatrix} >>> unmirror0(array([1, 2, 3, 2])) array([1, 2, 3]) ''' return f[:(len(f)/2+1)]
[docs]def unmirror1(f): r''' .. math:: \mathbf{M}_{1}^{\dagger}:\quad\begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-1}\\ x_{N}\\ \vdots\\ x_{2N-1} \end{bmatrix}\mapsto\begin{bmatrix}x_{0}\\ x_{1}\\ \vdots\\ x_{N-1} \end{bmatrix} >>> unmirror1(array([1, 2, 3, -3, -2, -1])) array([1, 2, 3]) ''' return f[:(len(f)/2)] ######### # Methods using FFT #########
def H_exp(a): r''' :math:`\int_{x-h/2}^{x+h/2}f(\xi) \exp(i \xi) d\xi` This transformation is singular non-invertible. ''' N = len(a) h = 2*pi/N c = I_diag(N+2, -h/2, h/2) a = F(a) b = roll(c[2:]*a,+1) return Finv(b) def H_sin(a): r''' :math:`\int_{x-h/2}^{x+h/2}f(\xi) \sin(\xi) d\xi` This transformation is singular non-invertible. >>> x = linspace(0, 2*pi, 6)[:-1] >>> round_(real( H_sin(sin(2*x)) ), 3) array([ 0.271, 0.438, -0.573, -0.573, 0.438]) ''' N = len(a) h = 2*pi/N a = F(a) c = I_diag(N+2, -h/2, h/2) b = (roll(c[2:]*a,+1) - roll(c[:-2]*a,-1))/2j return Finv(b) def H_sin_dual(a): r''' :math:`\int_{x-h/2}^{x+h/2}f(\xi) \sin(\xi+\frac{h}{2}) d\xi` This transformation is singular non-invertible. >>> x = linspace(0, 2*pi, 6)[:-1] >>> round_(real( H_sin_dual(sin(2*x)) ), 3) array([ 0.219, 0.573, -0.084, -0.844, 0.135]) ''' N = len(a) h = 2*pi/N a = F(a) c = I_diag(N+2, -h/2, h/2) s = exp(1j*h/2) b = (roll(c[2:]*a,+1)*s - roll(c[:-2]*a,-1)/s)/2j return Finv(b)
[docs]def Omega(N): r''' .. math:: \mathbf{\Omega}_{nn}=\sin\left(nh\right) ''' h = 2*pi/N o = sin( arange(N)*h ) return o
[docs]def Omega_d(N): r''' .. math:: \mathbf{\widetilde{\Omega}}_{nn} =\sin\left(\left(n+\frac{1}{2}\right)h\right) If :math:`\widetilde{\omega}` is the length of each dual edge, then .. math:: \mathbf{\widetilde{\Omega}} =\text{diag}(\mathbf{M}_{1}^{\dagger} \mathcal{F}^{-1}{\mathbf{I}^{-\frac{h}{2},\frac{h}{2}}}^{-1} \mathcal{F}\mathbf{M}_{1}^{-}\widetilde{\omega}) ''' h = 2*pi/N o = sin( (arange(N)+0.5)*h ) return o
[docs]def half_edge_base(N): r''' This is the discrete version of :math:`\delta(x)` - the basis functions for the half edge at -1 .. math:: \mathbf{B}= \text{diag}\left(\underset{N}{\underbrace{ (N-1)^{2},-\frac{1}{2},\frac{1}{2},-\frac{1}{2},\cdots} }\right) .. math:: \mathbf{B}^\dagger= \text{diag}\left(\underset{N}{\underbrace{ \cdots,-\frac{1}{2},\frac{1}{2},-\frac{1}{2},(N-1)^{2}} }\right) >>> half_edge_base(3) array([ 4.5, -0.5, 0.5]) ''' a = 0.5*ones(N) a[1::2] += -1 a[0] += (N-1)**2 return a
def half_edge_integrals(f): ''' >>> approx( half_edge_integrals([1,0,0]), gamma(3, 0)) True >>> approx( half_edge_integrals([0,1,0]), gamma(3, 1)) True ''' N = len(f) h = 2*pi/N f = Hinv(f) p = I(f, 0, h/2) return p[0]
[docs]def pick(f, n): r''' Pick the nth element in the array f. .. math:: \mathcal{\mathbf{P}}^{n}:\quad\begin{bmatrix}x_{0}\\ \vdots\\ x_{n-1}\\ x_{n}\\ x_{n+1}\\ \vdots\\ x_{N-1} \end{bmatrix}\mapsto\begin{bmatrix}x_{n}\\ \vdots\\ x_{n}\\ x_{n}\\ x_{n}\\ \vdots\\ x_{n} \end{bmatrix} .. math:: \mathcal{\mathbf{P}}^{0}=\begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 1\\ 1\\ 1\\ 1 \end{pmatrix},\quad\mathcal{\mathbf{P}}^{1}=\begin{pmatrix} 0 & 1 & 0 & 0 & 0\\ & 1\\ & 1\\ & 1\\ & 1 \end{pmatrix},\dots ''' return f[n]*ones(f.shape[0])
def reverse(f): r''' Reverse array. .. math:: \mathbf{R}:\quad\begin{bmatrix}x_{0}\\ \vdots\\ x_{n-1}\\ x_{n}\\ x_{n+1}\\ \vdots\\ x_{N-1} \end{bmatrix}\mapsto\begin{bmatrix}x_{N-1}\\ \vdots\\ x_{n+1}\\ x_{n}\\ x_{n-1}\\ \vdots\\ x_{0} \end{bmatrix} .. math:: \mathbf{R} = \begin{pmatrix} & & & & 1\\ & & & 1\\ & & 1\\ & 1\\ 1 \end{pmatrix} ''' return f[::-1] def I_space(a, b): return lambda f: Finv(F(f)*I_diag(f.shape[0], a, b)) def I_space_inv(a, b): return lambda f: Finv(F(f)/I_diag(f.shape[0], a, b)) def T_space(a): return lambda f: Finv(F(f)*S_diag(f.shape[0], a)) def T_space_inv(a): return lambda f: Finv(F(f)/S_diag(f.shape[0], a)) def E_space(n): return lambda f: Finv(extend(F(f), n)) def E_space_inv(n): return lambda f: Finv(unextend(F(f), n)) def matA(a): return concatenate(( [a[0]], a[1:-1:2]+a[2:-1:2] , [a[-1]] )) def matB(f): b = half_edge_base(f.shape[0]) return b*f[0] def matB1(f): b = half_edge_base(f.shape[0]) return b[::-1]*f[-1] def matC(f): return concatenate(( [0], f[1:-1], [0] ))
[docs]def H0d_cheb(f): r''' .. math:: \widetilde{\mathbf{H}}^0 = {\mathbf{M}_1}^{\dagger} \mathbf{I}^{-\frac{h}{2}, \frac{h}{2}} \widetilde{\mathbf{\Omega}} \mathbf{M}_1^+ ''' f = mirror1(f, +1) N = f.shape[0]; h = 2*pi/N f = f*Omega_d(N) f = I_space(-h/2, h/2)(f) f = unmirror1(f) return f
[docs]def H1_cheb(f): r''' .. math:: \mathbf{H}^1 = {\mathbf{M}_1}^{\dagger} \widetilde{\mathbf{\Omega}}^{-1} {\mathbf{I}^{-\frac{h}{2}, \frac{h}{2}}}^{-1} \mathbf{M}_1^- ''' f = mirror1(f, -1) N = f.shape[0]; h = 2*pi/N f = I_space_inv(-h/2, h/2)(f) f = f/Omega_d(N) f = unmirror1(f) return f
[docs]def H0_cheb(f): ''' Attempt to simplify the hodge-star. ''' N = f.shape[0]; h = pi/(N-1) f = mirror0(f, +1) f = E_space(N-1)(f) f = I_space(0, h/2)(f*Omega(f.shape[0])) f = unmirror1(f) f = matA(f) return real(f)
def H0_cheb_alternate(f): r''' .. math:: \mathbf{H}^0 = \mathbf{M}_0^{\dagger} (\mathcal{A}^{0} \mathbf{E}^{-1} \mathbf{I}^{-\frac{h}{2}, 0} + \mathcal{A}^{N-1} \mathbf{E}^{-1} \mathbf{I}^{0, +\frac{h}{2}}) \mathbf{\Omega} \mathbf{E}^{1} \mathbf{M}_0^{+} ''' N = f.shape[0] f = mirror0(f, +1) h = 2*pi/f.shape[0] f = E_space(1)(f) f = f*Omega(f.shape[0]) l, r = I_space(-h/2, 0)(f), I_space(0, +h/2)(f) l, r = E_space_inv(1)(l), E_space_inv(1)(r) l[0], r[N-1] = 0, 0 f = l+r f = unmirror0(f) #TODO: Is it possible to avoid discarding the imaginary part? f = real(f) return f
[docs]def H1d_cheb(f): r''' .. math:: \widetilde{\mathbf{H}}^{1} = \mathbf{M}_{0}^{\dagger}\left(\mathbf{T^{-\frac{h}{2}}}\mathbf{\Omega}^{-1}\mathbf{T}^{\frac{h}{2}}-\mathbf{B}\mathbf{I}^{0,\frac{h}{2}}-\mathbf{B}^{\dagger}\mathbf{I}^{-\frac{h}{2},0}\right)\mathbf{I}^{-\frac{h}{2},\frac{h}{2}}^{-1}\mathbf{M}_{0}^{-}\mathbf{C}+\mathbf{B}+\mathbf{B}^{\dagger} ''' N = f.shape[0]; h = pi/(N-1) B = half_edge_base(N) def endpoints(f): f0 = mirror0(matC(f), -1) aa = f - unmirror0(I_space(0, h/2)(I_space_inv(-h/2, h/2)(f0))) bb = f - unmirror0(I_space(-h/2, 0)(I_space_inv(-h/2, h/2)(f0))) return matB(aa) + matB1(bb) def midpoints(f): f = mirror0(matC(f), -1) # Shift function with S, Sinv to avoid division by zero at x=0, x=pi f = I_space_inv(-h/2, h/2)(f) f = T_space(+h/2)(f) f = f/Omega_d(f.shape[0]) f = T_space(-h/2)(f) f = unmirror0(f) return f return midpoints(f) + endpoints(f)
[docs]def laplacian(g): ''' Laplacian Operator ''' D, DD = g.derivative() H0, H1, H0d, H1d = g.hodge_star() L = lambda x: H1d(DD(H1(D(x)))) Ld = lambda x: H1(D(H1d(DD(x)))) return L, Ld
[docs]class Grid_1D_Chebyshev: def __init__(self, n, xmin=-1, xmax=+1): N = n dimension = 1 # 2n-1 points: n primal, n-1 dual x = sin(linspace(-pi/2, pi/2, 2*n-1)) p = 0.5*(xmin*(1-x) + xmax*(1+x)) verts = p[::2] delta = diff(verts) edges = (verts[:-1], verts[1:]) verts_dual = p[1::2] tmp = concatenate(([p[0]], verts_dual, [p[-1]])) delta_dual = diff(tmp) edges_dual = (tmp[:-1], tmp[1:]) V = verts S0 = arange(len(V)) S1 = (S0[:-1], S0[1:]) Vd = verts_dual S0d = arange(len(Vd)) S1d = (S0d[:-1], S0d[1:]) self.dimension = dimension self.n = n self.xmin = xmin self.xmax = xmax self.delta = delta self.delta_dual = delta_dual self.pnts = p self.verts = verts self.verts_dual = verts_dual self.edges = edges self.edges_dual = edges_dual def projection(self): P0 = lambda f: f(self.verts) P1 = lambda f: integrate_chebyshev(self.verts, f) P0d = lambda f: f(self.verts_dual) P1d = lambda f: integrate_chebyshev_dual( concatenate(([-1], self.verts_dual, [+1])), f) return P0, P1, P0d, P1d def basis_fn(self): n = self.n B0 = [lambda x, i=i: psi0(n, i, x) for i in range(n)] B1 = [lambda x, i=i: psi1(n, i, x) for i in range(n-1)] B0d = [lambda x, i=i: psid0(n, i, x) for i in range(n-1)] B1d = [lambda x, i=i: psid1(n, i, x) for i in range(n)] return B0, B1, B0d, B1d def reconstruction(self): R0, R1, R0d, R1d = reconstruction(self.basis_fn()) return R0, R1, R0d, R1d def boundary_condition(self, f): bc = zeros((self.n, )) bc[ 0] = -f(self.xmin) bc[-1] = +f(self.xmax) return bc def derivative(self): D0 = lambda f: diff(f) D0d = lambda f: diff(concatenate(([0], f, [0]))) return D0, D0d def hodge_star(self): H0 = H0_cheb H1 = H1_cheb H0d = H0d_cheb H1d = H1d_cheb return H0, H1, H0d, H1d def contraction(self, V): ''' Implement contraction where V is the one-form corresponding to the vector field. ''' def C1(alpha): return None return C1