Source code for geotecha.mathematics.hankel

# geotecha - A software suite for geotechncial engineering
# Copyright (C) 2018  Rohan T. Walker (rtrwalker@gmail.com)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see http://www.gnu.org/licenses/gpl.html.
"""Hankel transforms."""


from __future__ import division, print_function
import matplotlib.pyplot
import numpy as np
from scipy import integrate
from scipy.special import jn_zeros
from scipy.special import jn
from matplotlib import pyplot as plt
import functools

import unittest
from numpy.testing import assert_allclose
from numpy.polynomial.polynomial import Polynomial

from geotecha.mathematics.quadrature import gl_quad
from geotecha.mathematics.quadrature import gk_quad
from geotecha.mathematics.quadrature import gauss_kronrod_abscissae_and_weights
from geotecha.mathematics.quadrature import gauss_legendre_abscissae_and_weights
from geotecha.mathematics.quadrature import shanks_table
from geotecha.mathematics.quadrature import shanks

[docs]class HankelTransform(object): """Hankel transform of integer order Applies Gauss-Kronrod quadrature between zeros of the integrand. Then uses shanks extrapolation. Parameters ---------- func : function or method Function to apply hankel trasnform to. `func` is called with func(r, *args). `func` should be vectorised in r. args : tuple, optional Arguments to pass to f. Default args=(). order : integer, optional Order of Hankel transform. Default order=0. m : int, optional Number of segments to break the integration interval into. Each segment will be between the zeros of the bessel function or order= `order' divide by r Default ,=20. points : list/array of float, optional Points in addition to those defined by m at which to split the integral. Use this if `f` itself oscillates or there are discontinuities. Points will only be included that are less than the `m` zeros mentined above. Default points=None i.e. no extra points. Basically the function is never evaluated at any ofthe points, rather they form the boundary four gauss quadrature. ng : [7,10,15,20,25,30], optional Number of gauss points to use in integration after first zero. Default ng=10. Number of Kronrod points will automatically by 2 * ng + 1. ng0 : [7,10,15,20,25,30], optional Number of gauss points to use in integrating between 0 and first zero. Default ng0=20. shanks_ind : int, optional Start position of intervals (not including the first interval) from which to begin shanks transformation. Default shanks_ind=None i.e. no extrapolation. The first interval will never be included. Be careful when using shanks extrapolation; make sure you only begin to use it after the intgrand is well behaved. Use the plot_integrand method to check your integrand. Returns ------- f : float Value of transform at r. Notes ----- The Hankel Transform of order :math:`\\nu` is given by: .. math:: F_\\nu(s)=\mathcal{H}_\\nu\\{f(r)\\} = \\int_0^{\\infty}rf(r)J_\\nu(sr)\\,\\mathrm{d}r Provided :math:`\\nu\\gt1/2` the inverse hankel transform is the same as the normal transform: .. math:: f(r)=\mathcal{H}_{\\nu}^{-1}\\{F_{\\nu}(s)\\} = \\int_0^{\\infty}sF_\\nu(s)J_\\nu(sr)\\,\\mathrm{d}s References ---------- .. [1] Piessens, Robert. 2000. 'Chapter 9 The Hankel Transform'. In The Transforms and Applications Handbook, edited by Alexander D. Poularikas, 2nd edition. Boca Raton, USA: CRC Press. """ def __init__(self, func, args=(), order=0, m=20, points=None, ng=10, ng0=20, shanks_ind=None): self.func = func self.args = args self.order = order self.m = m if points is None: self.points = None else: self.points = np.atleast_1d(points) self.ng = ng self.ng0 = ng0 self.shanks_ind = shanks_ind self.zeros_of_jn()
[docs] def zeros_of_jn(self): """Roots of Jn for determining integration intervals (0 prepended)""" self.jn_0s = np.zeros(self.m + 1, dtype=float) self.jn_0s[1:] = jn_zeros(self.order, self.m) return
def _integrand(self, s, r, *args): """ Parameters ---------- s : float transform variable. r : float independent variable. i.e. integrate w.r.t. to r """ return self.func(r, *args) * r * jn(self.order, s * r) def __call__(self, s, a=0, b=np.inf): """transform f(r, *args) at s Parameters ---------- s : scalar transform variable, i.e. coordinate to evaluate transform at a, b : float, optional limits of integration. default a=0, b=np.inf. A hankel transform is usually from 0 to infinity. However, if you have a function that is defined on [a,b] and zero elsewhere then you can restrict integration to that interval (shanks_ind will be ignored if b!=np.inf). Returns ------- F : float Transformed functin evaluated at s. err_est : float Error estimate. For each interval (i.e. between bessel zeros and any specified points) sum up 200*abs(G-K)**1.5. The error is calculated before any shanks extrapolation so the error is just a measure of the difference between the coarse Gauss quadrature and the finer Kronrod quadrature. """ integrand = functools.partial(self._integrand, s) zeros = self.jn_0s / s if not self.points is None: #zeros becomes both zeros and interesting points/singularities zeros = np.unique(list(zeros) + list(self.points[(self.points < zeros[-1]) & (self.points>0)])) if (a!=0) or (b!=np.inf): zeros = np.unique(zeros.clip(a, b)) #1st segment igral0, err_est0 = gk_quad(integrand, zeros[0], zeros[1], self.args, self.ng0) #remaining segments if len(zeros)>2: igralm, err_estm = gk_quad(integrand, zeros[1:-1], zeros[2:], self.args, self.ng) else: return igral0[0], err_est0[0] if (self.shanks_ind is None) or (b!=np.inf): igral = igral0 + np.sum(igralm) else: igralm.cumsum(out=igralm) igral = igral0 + shanks(igralm, self.shanks_ind) err_est = (200*np.abs(err_est0))**1.5 + np.sum((200*np.abs(err_estm))**1.5) return igral[0], err_est[0]
[docs] def plot_integrand(self, s, npts = 1000): """Plot the integrand Parameters ---------- s : float Transform variable, i.e. point to evaluate transform at npts : int, optional Number of points to plot. Default npts=1000. Returns ------- fig : matplotlib.Figure Use plt.show to plot Notes ----- Use this to check if your parameters are appropriate. """ integrand = functools.partial(self._integrand, s) zeros=self.jn_0s / s if not self.points is None: zeros = np.unique(list(zeros) + list(self.points[(self.points < zeros[-1]) & (self.points>0)])) x = np.linspace(0.001, zeros[-1], npts) # y = self.f(x, *self.args) * x * jn(self.order, x * r) y = integrand(x, *self.args) fig = plt.figure() ax = fig.add_subplot('111') ax.plot(x,y, '-') ax.plot(zeros, np.zeros_like(zeros), 'ro', ms=5) # ax.plot(x, x*jn(self.order, r*x), label='rJn(r*s)') ax.set_title('$rf(r)J_\\nu(sr)\\,dr)$ at s={:.3g}, f(r) is {}(r, *{})'.format(s, self.func.__name__, self.args)) ax.set_xlabel('x') ax.set_ylabel('y') ax.grid() return fig
[docs]def vhankel_transform(f, r, args=(), order=0, m=20, ng=20, shanks_ind=None): """Hankel transform of f(r) This is a vectorised Hankel transform Parameters ---------- f : function or method Function to apply hankel trasnform to. f is called with f(s, *args). r : 1d array Coordinate(s) to evaluate transform at. args : tuple, optional Arguments to pass to f, default args=(). order : integer, optional Order of hankel transform. Default order=0. m : int, optional Number of segments to break the integration interval into. Each segment will be between the zeros of the bessel function. Default m=20. ng : [2-20, 32, 64, 100], optional Number of gauss points to use in integration. Default ng=20. shanks_ind : int, optional Start position of intervals to start shanks extrapolation. Default shanks_ind=None i.e. no extrapolation. Be careful when using shanks extrapolation; make sure you only begin to use it after the intgrand is well behaved. Returns ------- f : 1d array of float Value of transform at r. Notes ----- The Hankel Transform of order :math:`\\nu` is given by: .. math:: F_\\nu(s)=\mathcal{H}_\\nu\\{f(r)\\} = \\int_0^{\\infty}rf(r)J_\\nu(sr)\\,\\mathrm{d}r Provided :math:`\\nu\\gt1/2` the inverse hankel transform is the same as the normal transform: .. math:: f(r)=\mathcal{H}_{\\nu}^{-1}\\{F_{\\nu}(s)\\} = \\int_0^{\\infty}sF_\\nu(s)J_\\nu(sr)\\,\\mathrm{d}s Note that because this implementation does not allow for input of extra point to break up the integration inteval, there is no way to account for singularities and other oscillations. If you need this control then see the HankelTransform class which is not vectorized but provides a few more options. See Also -------- HankelTransform : Non vectorised Hankel transform. References ---------- .. [1] Piessens, Robert. 2000. 'Chapter 9 The Hankel Transform'. In The Transforms and Applications Handbook, edited by Alexander D. Poularikas, 2nd edition. Boca Raton, USA: CRC Press. """ ri = np.atleast_1d(r) xk_, wk = gauss_legendre_abscissae_and_weights(ng) # integration intervals zeros = np.zeros(m + 1, dtype=float) zeros[1:] = jn_zeros(order, m) aj = zeros[0:-1] bj = zeros[1:] ri = ri[:, np.newaxis, np.newaxis] aj = aj[np.newaxis, :, np.newaxis] bj = bj[np.newaxis, :, np.newaxis] xk_ = xk_[np.newaxis, np.newaxis, :] wk = wk[np.newaxis, np.newaxis, :] aij = aj / ri bij = bj / ri bma = (bij - aij) / 2 # b minus a bpa = (aij + bij) /2 # b plus a xijk = bma * xk_ + bpa # xj_ are in [-1, 1] so need to transform to [a, b] fijk = f(xijk, *args) fijk *= jn(order, ri * xijk) fijk *= xijk igral = bma[:,:,0] * np.sum(fijk * wk, axis=2) if shanks_ind is None: return igral.sum(axis=1) else: #extrapolate igral.cumsum(axis=1 , out=igral) return shanks(igral, shanks_ind)
##Hankel transform pairs ##zero order #def hankel1(s, a): # """a/(s**2 + a**2)**1.5""" # #H(hankel1)=exp(-a* r) # return a/(s**2 + a**2)**1.5 #def hankel1_(r, a): # """exp(-a*r)""" # return np.exp(-a*r) # #def hankel2(s, *args): # "1/s" # #H(hankel2)=1/r # return 1/s #def hankel2_(r, *args): # "1/r" # return 1/r # #def hankel3(s,a): # "1/s*jn(0,a/s)" # #H(hankel3)=1/rJ0(2*(a*r)**0.5) # return 1/s*jn(0,a/s) #def hankel3_(r, a): # "1/rJ0(2*(a*r)**0.5)" # return 1/r*jn(0,(2*(a*r)**0.5)) # ##integer order #def hankel4(s, a, v=0): # """(sqrt(s**2+a**2)-a)**v/(s**v*sqrt(s**2+a**2))""" # #H(hankel4)=exp(-a*r)/r # return (np.sqrt(s**2 + a**2) - a)**v/(s**v*np.sqrt(s**2+a**2)) # #def hankel4_(r, a, *args): # """exp(-a*r)/r""" # return np.exp(-a*r)/r # #def hankel5(s, a, v=0): # """s**v/(2*a**2)**(v+1)*exp(-s**2/(4*a**2))""" # #H(hankel5)=exp(-a**2*r**2)*r**v # return s**v/(2*a**2)**(v+1)*np.exp(-s**2/(4*a**2)) #def hankel5_(r, a, v=0): # """exp(-a**2*r**2)*r**v""" # return np.exp(-a**2*r**2)*r**v
[docs]def scratch(): """ """
if __name__ == '__main__': import nose nose.runmodule(argv=['nose', '--verbosity=3', '--with-doctest']) # nose.runmodule(argv=['nose', '--verbosity=3'])