Merge pull request #6 from prmiles/calc_fd

Added FD calculator

Author: Paul Miles

pyfod/GaussLaguerre.py

@@ -7,16 +7,20 @@
 """
 
 import numpy as np
-import sys
 from sympy.integrals.quadrature import gauss_gen_laguerre as sp_gauss_laguerre
 import sympy as sp
+from pyfod.utilities import check_function
 
 
 class GaussLaguerre:
 
     def __init__(self, N=5, start=0.0, finish=1.0, alpha=0.0,
                  f=None, extend_precision=True, n_digits=30):
         self.description = 'Gaussian-Laguerre Quadrature'
+        self.start = start
+        self.finish = finish
+        self.alpha = alpha
+        self.f = f
         if extend_precision is False:
             points, weights = np.polynomial.laguerre.laggauss(deg=N)
             self.points = 1 - np.exp(-points)
@@ -29,73 +33,100 @@ def __init__(self, N=5, start=0.0, finish=1.0, alpha=0.0,
                     np.ones(shape=len(points))) - points.applyfunc(sp.exp)
             weights = sp.Array(weights)
         self.weights = weights
-        self.start = start
-        self.finish = finish
-        self.alpha = alpha
-        self.f = f
+        self.initial_weights = weights.copy()
+        self.update_weights(alpha=alpha)
 
-    def integrate(self, f=None, alpha=None):
-        if alpha is None:
-            alpha = self.alpha
-        if f is None:
-            f = self.f
-        if f is None:
-            sys.exit('No function defined... provide function f')
+    def integrate(self, f=None):
+        f = check_function(f, self.f)
+        self.f = f
         # transform kernel
         span = self.finish - self.start
         # check if sympy
         if isinstance(self.points, sp.Array):
             evalpoints = self.points.applyfunc(
                     lambda x: span*x + self.start)
             feval = evalpoints.applyfunc(f)
-            coef = self.points.applyfunc(
-                    lambda x: span**(1-alpha)*(1-x)**(-alpha))
             s = 0
-            for ii, (w, f, t) in enumerate(zip(self.weights, feval, coef)):
-                s += w*f*t
+            for ii, (w, f) in enumerate(zip(self.weights, feval)):
+                s += w*f
             return s
         else:
-            coef = span**(1-alpha)*(1-self.points)**(-alpha)
-            s = (self.weights*(coef*f(span*self.points + self.start))).sum()
+            s = (self.weights*(f(span*self.points + self.start))).sum()
             return s
 
+    def update_weights(self, alpha=None):
+        if alpha is None:
+            alpha = self.alpha
+        self.alpha = alpha
+        span = self.finish - self.start
+        # check if sympy
+        if isinstance(self.points, sp.Array):
+            coef = self.points.applyfunc(
+                    lambda x: span**(1-alpha)*(1-x)**(-alpha))
+            wtmp = []
+            for ii, (c, w) in enumerate(zip(coef, self.initial_weights)):
+                wtmp.append(c*w)
+            self.weights = sp.Array(wtmp)
+        else:
+            coef = span**(1-alpha)*(1-self.points)**(-alpha)
+            self.weights = self.initial_weights*coef
+
 
 if __name__ == '__main__':  # pragma: no cover
 
     n_digits = 50
-    alpha = 0.9
     start = 0.0
     finish = 1.0
-    dt = 1e-6
     N = 24
+    '''
+    Test normal precision
+    '''
+    # f(t) = exp(2t)
+
+    def fexpnp(x):
+        return np.exp(2*x)
 
     def fcosnp(x):
-        return np.cos(2*x) + 3
+        return np.cos(2*x)
 
+    # Test alpha = 0.0
+    alpha = 0.0
     GLag = GaussLaguerre(N=N, start=start, finish=finish,
-                         f=fcosnp, alpha=0, extend_precision=False)
+                         f=fexpnp, alpha=alpha, extend_precision=False)
     F1 = GLag.integrate()
-    print('Int(cos(2t)+3) = {} ({})'.format(F1, 3.4546))
+    F2 = GLag.integrate(f=fcosnp)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.19453))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.454649))
+
+    '''
+    Test extended precision
+    '''
 
     def fexp(x):
         return sp.exp(2*x)
-    GLag = GaussLaguerre(N=N, start=start, finish=finish,
-                         f=fexp, alpha=alpha, n_digits=n_digits)
-    F1 = GLag.integrate()
-    GLag = GaussLaguerre(N=N, start=start, finish=finish-dt,
-                         f=fexp, alpha=alpha, n_digits=n_digits)
-    F2 = GLag.integrate()
-    print('D[f(t)=exp(2t)] = {} ({})'.format(
-            (F1-F2)/(dt*sp.gamma(1-0.9)), 13.815))
 
     def fcos(x):
         return sp.cos(2*x)
 
+    # Test alpha = 0.0
+    alpha = 0.0
     GLag = GaussLaguerre(N=N, start=start, finish=finish,
-                         f=fcos, alpha=alpha, n_digits=n_digits)
-    F1 = GLag.integrate()
-    GLag = GaussLaguerre(N=N, start=start, finish=finish-dt,
-                         f=fcos, alpha=alpha, n_digits=n_digits)
-    F2 = GLag.integrate()
-    print('D[f(t)=cos(2t)] = {} ({})'.format(
-            (F1-F2)/(dt*sp.gamma(1-0.9)), -1.779))
+                         alpha=alpha, n_digits=n_digits)
+    F1 = GLag.integrate(f=fexp)
+    F2 = GLag.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.19453))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.454649))
+    # Test alpha = 0.1
+    alpha = 0.1
+    GLag.update_weights(alpha=alpha)
+    F1 = GLag.integrate(f=fexp)
+    F2 = GLag.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.749))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.457653))
+    # Test alpha = 0.9
+    alpha = 0.9
+    GLag.update_weights(alpha=alpha)
+    F1 = GLag.integrate(f=fexp)
+    F2 = GLag.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 65.2162))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, -2.52045))

pyfod/GaussLegendre.py

@@ -7,21 +7,41 @@
 """
 
 import numpy as np
+from pyfod.utilities import check_alpha
+from pyfod.utilities import check_function, check_singularity
 
 
 class GaussLegendre:
 
-    def __init__(self, N=5, start=0.0, finish=1.0):
+    def __init__(self, N=5, start=0.0, finish=1.0,
+                 alpha=0.0, f=None, singularity=None):
         self.description = 'Gaussian-Legendre Quadrature'
+        check_alpha(alpha)
         h = (finish - start)/N
-        self.points = self.gauss_points(N=N, h=h, start=start)
-        self.weights = self.gauss_weights(N=N, h=h)
-
-    def integrate(self, f):
+        self.alpha = alpha
+        self.f = f
+        self.finish = finish
+        self.singularity = check_singularity(singularity, self.finish)
+        self.points = self._gauss_points(N=N, h=h, start=start)
+        self.weights = self._gauss_weights(N=N, h=h)
+        self.initial_weights = self.weights.copy()
+        self.update_weights(alpha=alpha)
+
+    def update_weights(self, alpha=None):
+        if alpha is None:
+            alpha = self.alpha
+        self.alpha = alpha
+        # update weights based on alpha
+        self.weights = self.initial_weights*(
+                self.singularity - self.points)**(-alpha)
+
+    def integrate(self, f=None):
+        f = check_function(f, self.f)
+        self.f = f
         return (self.weights*f(self.points)).sum()
 
     @classmethod
-    def base_gauss_points(cls):
+    def _base_gauss_points(cls):
 
         # base points
         gpts = np.zeros([4])
@@ -32,7 +52,7 @@ def base_gauss_points(cls):
         return gpts
 
     @classmethod
-    def base_gauss_weights(cls, h):
+    def _base_gauss_weights(cls, h):
         # define the Gauss weights for a four point quadrature rule
         w = np.zeros([4])
         w[0] = 49*h/(12*(18 + np.sqrt(30)))
@@ -42,7 +62,7 @@ def base_gauss_weights(cls, h):
         return w
 
     @classmethod
-    def interval_gauss_points(cls, base_gpts, N, h, start):
+    def _interval_gauss_points(cls, base_gpts, N, h, start):
         # determines the Gauss points for all N intervals.
         Gpts = np.zeros([4*N])
         for gct in range(N):
@@ -51,16 +71,16 @@ def interval_gauss_points(cls, base_gpts, N, h, start):
                                        + base_gpts[ell]*h + start)
         return Gpts
 
-    def gauss_points(self, N, h, start):
+    def _gauss_points(self, N, h, start):
         # base points
-        gpts = self.base_gauss_points()
+        gpts = self._base_gauss_points()
         # determines the Gauss points for all N intervals.
-        Gpts = self.interval_gauss_points(gpts, N, h, start)
+        Gpts = self._interval_gauss_points(gpts, N, h, start)
         return Gpts
 
-    def gauss_weights(self, N, h):
+    def _gauss_weights(self, N, h):
         # determine the Gauss weights for a four point quadrature rule
-        w = self.base_gauss_weights(h)
+        w = self._base_gauss_weights(h)
         # copy the weights to form a vector for all N intervals
         weights = w.copy()
         for gct in range(N-1):
@@ -70,10 +90,35 @@ def gauss_weights(self, N, h):
 
 if __name__ == '__main__':  # pragma: no cover
 
-    GQ = GaussLegendre(N=10, start=1.0, finish=12.0)
-
-    def f(t):
-        return np.cos(2*t) + 3
-
-    a = GQ.integrate(f)
-    print('a = {}'.format(a))
+    start = 0.0
+    finish = 1.0
+    N = 2400
+    GQ = GaussLegendre(N=N, start=start, finish=finish)
+
+    def fexp(x):
+        return np.exp(2*x)
+
+    def fcos(x):
+        return np.cos(2*x)
+
+    # Test alpha = 0.0
+    alpha = 0.0
+    GLeg = GaussLegendre(N=N, start=start, finish=finish, alpha=alpha)
+    F1 = GLeg.integrate(f=fexp)
+    F2 = GLeg.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.19453))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.454649))
+    # Test alpha = 0.1
+    alpha = 0.1
+    GLeg.update_weights(alpha=alpha)
+    F1 = GLeg.integrate(f=fexp)
+    F2 = GLeg.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.749))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.457653))
+    # Test alpha = 0.9
+    alpha = 0.9
+    GLeg.update_weights(alpha=alpha)
+    F1 = GLeg.integrate(f=fexp)
+    F2 = GLeg.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 65.2162))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, -2.52045))

pyfod/GaussLegendreRiemannSum.py

@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 25 14:51:32 2019
+
+@author: prmiles
+"""
+
+from pyfod.GaussLegendre import GaussLegendre
+from pyfod.RiemannSum import RiemannSum
+from pyfod.utilities import check_function
+import numpy as np
+
+
+class GaussLegendreRiemannSum(object):
+
+    def __init__(self, NGQ=5, NRS=20, pGQ=0.9,
+                 start=0.0, finish=1.0, alpha=0.0, f=None):
+        self.description = 'Gaussian Quadrature, Riemann-Sum'
+        # setup GQ points/weights
+        switch_time = (finish - start)*pGQ
+        self.GQ = GaussLegendre(N=NGQ, start=start, finish=switch_time,
+                                alpha=alpha, singularity=finish, f=f)
+        # setup RS points/weights
+        self.RS = RiemannSum(N=NRS, start=switch_time,
+                             finish=finish, alpha=alpha, f=f)
+        self.alpha = alpha
+        self.pGQ = pGQ
+        self.f = f
+
+    def integrate(self, f=None):
+        f = check_function(f, self.f)
+        self.f = f
+        return self.GQ.integrate(f=f) + self.RS.integrate(f=f)
+
+    def update_weights(self, alpha=None):
+        if alpha is None:
+            alpha = self.alpha
+        self.alpha = alpha
+        self.GQ.update_weights(alpha=alpha)
+        self.RS.update_weights(alpha=alpha)
+
+
+if __name__ == '__main__':  # pragma: no cover
+
+    start = 0.0
+    finish = 1.0
+    NGQ = 100
+    NRS = 2400
+    pGQ = 0.95
+
+    def fexp(x):
+        return np.exp(2*x)
+
+    def fcos(x):
+        return np.cos(2*x)
+
+    # Test alpha = 0.0
+    alpha = 0.0
+    GRS = GaussLegendreRiemannSum(pGQ=pGQ, NGQ=NGQ, NRS=NRS, start=start,
+                                  finish=finish)
+    F1 = GRS.integrate(f=fexp)
+    F2 = GRS.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.19453))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.454649))
+    # Test alpha = 0.1
+    alpha = 0.1
+    GRS.update_weights(alpha=alpha)
+    F1 = GRS.integrate(f=fexp)
+    F2 = GRS.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 3.749))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, 0.457653))
+    # Test alpha = 0.9
+    alpha = 0.9
+    GRS.update_weights(alpha=alpha)
+    F1 = GRS.integrate(f=fexp)
+    F2 = GRS.integrate(f=fcos)
+    print('Int(exp(2t)/(1-t)^{}) = {} ({})'.format(alpha, F1, 65.2162))
+    print('Int(cos(2t)/(1-t)^{}) = {} ({})'.format(alpha, F2, -2.52045))