Tweak approximation methods for consistency

lcapy/approximate.py

@@ -1,24 +1,20 @@
-"""This module contains functions for approximating expressions.
+"""This module contains functions for approximating SymPy expressions.
 
-Copyright 2021--23 Michael Hayes, UCECE
+Copyright 2021--2024 Michael Hayes, UCECE
 
 """
 
 from .simplify import expand_hyperbolic_trig
-from sympy import exp, diff
+from sympy import exp, diff, eye, zeros, Matrix, Poly
 from math import factorial
 
 
-def approximate_fractional_power(expr, method='pade', order=2):
+def approximate_fractional_power(expr, var, ndegree=2, ddegree=2):
     """This is an experimental method to approximate s**a, where a is
-    fractional, with a rational function using a Pade approximant.
+    fractional, with a rational function using a Pade approximant
+    of order [ndegree, ddegree]."""
 
-    """
-
-    if method != 'pade':
-        raise ValueError('Method %s unsupported, must be pade')
-
-    v = expr.var
+    v = var
 
     def query(expr):
 
@@ -51,39 +47,44 @@ def value2(expr):
         d = v**2 * (a**2 - 3 * a + 2) + v * (8 - a**2) + (a**2 + 3 * a + 2)
         return n / d
 
-    if order == 1:
+    if ndegree != ddegree:
+        raise ValueError('Require ndegree == ddegree')
+
+    if ndegree == 1:
         value = value1
-    elif order == 2:
+    elif ndegree == 2:
         value = value2
     else:
-        raise ValueError('Can only handle order 1 and 2 at the moment')
+        raise ValueError('Can only handle degree 1 and 2 at the moment')
 
-    expr = expr.expr
     expr = expr.replace(query, value)
 
     return expr
 
 
-def _approximate_exp_pade(expr, order=1, numer_order=None):
+def approximate_exp(expr, ndegree=1, ddegree=1):
+    """Approximate exp(a) with Pade approximant.  The best time-domain
+    response (without a jump) is achieved with 'ndegree == ddegree -
+    1'.  The best frequency-domain response is achieved with ndegree
+    == ddegree.
 
-    if numer_order is None:
-        numer_order = order
+    """
 
     def query(expr):
         return expr.is_Function and expr.func == exp
 
     def value(expr):
         arg = expr.args[0]
 
-        if order == 1 and numer_order == 1:
+        if ddegree == 1 and ndegree == 1:
             return (2 + arg) / (2 - arg)
-        elif order == 2 and numer_order == 2:
+        elif ddegree == 2 and ndegree == 2:
             return (12 + 6 * arg + arg**2) / (12 - 6 * arg + arg**2)
 
         from math import factorial
 
-        m = numer_order
-        n = order
+        m = ndegree
+        n = ddegree
         numer = 0
         denom = 0
 
@@ -102,24 +103,11 @@ def value(expr):
     return expr.replace(query, value)
 
 
-def approximate_exp(expr, method='pade', order=1, numer_order=None):
-    """Approximate exp(a).  The best time-domain response (without a jump)
-    is achieved with 'numer_order == order - 1'.  The best
-    frequency-domain response is achieved with numer_order ==
-    order."""
-
-    if method != 'pade':
-        raise ValueError('Method %s unsupported, must be pade')
-
-    return _approximate_exp_pade(expr, order, numer_order)
-
-
-def approximate_hyperbolic_trig(expr, method='pade', order=1, numer_order=None):
+def approximate_hyperbolic_trig(expr, ndegree=1, ddegree=1):
     """Approximate cosh(a), sinh(a), tanh(a)."""
 
     expr = expand_hyperbolic_trig(expr)
-    return approximate_exp(expr, method=method, order=order,
-                           numer_order=numer_order)
+    return approximate_exp(expr, ndegree=ndegree, ddegree=ddegree)
 
 
 def approximate_dominant_terms(expr, defs, threshold=0.01):
@@ -165,28 +153,83 @@ def value(expr):
     return expr.replace(query, value)
 
 
-def approximate_order(expr, var, order):
-    """Approximate expression by reducing order of polynomial to
-    specified order."""
+def approximate_degree(expr, var, degree):
+    """Approximate expression by reducing degree of polynomial to
+    specified degree."""
 
     from sympy import O
 
-    return (expr.expand() + O(var ** (order + 1))).removeO()
+    return (expr.expand() + O(var ** (degree + 1))).removeO()
+
+
+def approximate_taylor_coeffs(expr, var, degree=1, var0=0):
+
+    coeffs = [expr.subs(var, var0)]
+
+    prev = expr
+
+    for n in range(1, degree + 1):
+
+        prev = diff(prev, var)
+
+        coeffs.append(prev.subs(var, var0) / factorial(n))
+
+    return coeffs
 
 
-def approximate_taylor(expr, var, var0, degree):
+def approximate_taylor(expr, var, degree=1, var0=0):
     """Approximate expression using a Taylor series
     around `var = var0` to degree `degree`."""
 
     result = expr.subs(var, var0)
 
     prev = expr
 
-    for n in range(degree):
+    for n in range(1, degree + 1):
 
         prev = diff(prev, var)
 
         result += prev.subs(var, var0) * \
-            (var - var0)**(n + 1) / factorial(n + 1)
+            (var - var0)**n / factorial(n)
 
     return result
+
+
+def approximate_pade_coeffs_from_taylor_coeffs(a, m):
+
+    N = len(a) - 1
+    n = N - m
+    if n < 0:
+        raise ValueError('Degree %d must be smaller than %d' % (m, N))
+    A = eye(N + 1, n + 1)
+    B = zeros(N + 1, m)
+
+    for col in range(m):
+        for row in range(col + 1, N + 1):
+            B[row, col] = -a[row - col - 1]
+
+    C = A.hstack(A, B)
+    pq = C.solve(Matrix(a))
+    p = pq[:n + 1]
+
+    q = [1]
+    for col in range(1, m + 1):
+        q.append(pq[n + col])
+
+    return p[::-1], q[::-1]
+
+
+def approximate_pade_coeffs(expr, var, ndegree=2, ddegree=2):
+    """Approximate expression using a Pade series with numerator degree
+    `ndegree` and denominator degree `ddegree`."""
+
+    coeffs = approximate_taylor_coeffs(expr, var, ndegree + ddegree)
+
+    return approximate_pade_coeffs_from_taylor_coeffs(coeffs, ddegree)
+
+
+def approximate_pade(expr, var, ndegree=2, ddegree=2):
+
+    p, q = approximate_pade_coeffs(expr, var, ndegree, ddegree)
+
+    return Poly.from_list(p, gens=var) / Poly.from_list(q, gens=var)

lcapy/expr.py

@@ -40,7 +40,7 @@
 from .simplify import expand_hyperbolic_trig
 from .approximate import (approximate_fractional_power, approximate_exp,
                           approximate_hyperbolic_trig, approximate_dominant,
-                          approximate_order, approximate_taylor)
+                          approximate_degree, approximate_taylor, approximate_pade)
 from .config import heaviside_zero, unitstep_zero
 from collections import OrderedDict
 from warnings import warn
@@ -3558,97 +3558,122 @@ def approximate_dominant(self, defs, threshold=0.01):
         result = approximate_dominant(self.sympy, ndefs, threshold)
         return self.__class__(result, **self.assumptions)
 
-    def approximate_exp(self, method='pade', order=1, numer_order=None):
-        """Approximate exp(a).  The best time-domain response (without a jump)
-        is achieved with 'numer_order == order - 1', where `numer` is
-        the order of the numerator.  The best frequency-domain
-        response is achieved with `numer_order == order`."""
+    def approximate_exp(self, ndegree=2, ddegree=2):
+        """Approximate exp(a) with a Pade approximant order [ndegree,
+        ddegree].
 
-        expr = approximate_exp(self, method, order, numer_order)
+        The best time-domain response (without a jump) is achieved
+        with 'ndegree == ddegree - 1', where `numer` is the degree of
+        the numerator.  The best frequency-domain response is achieved
+        with `ndegree == ddegree`.
+
+        """
+
+        expr = approximate_exp(self, ndegree, ddegree)
         return self.__class__(expr, **self.assumptions)
 
-    def approximate_fractional_power(self, method='pade', order=2):
-        """This is an experimental method to approximate
-        s**a, where a is fractional, with a rational function using
-        a Pade approximant of order `order`."""
+    def approximate_fractional_power(self, ndegree=2, ddegree=2):
+        """Approximate s**a, where a is fractional, with a Pade approximant of
+        order [ndegree, ddegree].
 
-        expr = approximate_fractional_power(self, method, order)
+        """
+
+        expr = approximate_fractional_power(self.sympy, self.var,
+                                            ndegree, ddegree)
         return self.__class__(expr, **self.assumptions)
 
-    def approximate_hyperbolic_trig(self, method='pade', order=1,
-                                    numer_order=None):
-        """Approximate cosh(a), sinh(a), tanh(a) with a rational function using
-        a Pade approximant of order `order`."""
+    def approximate_hyperbolic_trig(self, ndegree=2, ddegree=2):
+        """Approximate cosh(a), sinh(a), tanh(a) with a Pade approximant of
+        order [ndegree, ddegree]."""
 
-        expr = approximate_hyperbolic_trig(self, method, order, numer_order)
+        expr = approximate_hyperbolic_trig(self, ndegree, ddegree)
         return self.__class__(expr, **self.assumptions)
 
-    def approximate_numer_order(self, order):
-        """Reduce order of numerator to `order`.
+    def approximate_numer_degree(self, degree):
+        """Reduce degree of numerator to `degree`.
         Warning, this assumes `self.var` is much smaller than 1."""
 
         N, D = self.as_N_D()
 
-        expr = approximate_order(N.sympy, self.var, order) / D.sympy
+        expr = approximate_degree(N.sympy, self.var, degree) / D.sympy
         return self.__class__(expr, **self.assumptions)
 
-    def approximate_denom_order(self, order):
-        """Reduce order of denominator to `order`.
+    def approximate_denom_degree(self, degree):
+        """Reduce degree of denominator to `degree`.
         Warning, this assumes `self.var` is much smaller than 1."""
 
         N, D = self.as_N_D()
 
-        expr = N.sympy / approximate_order(D.sympy, self.var, order)
+        expr = N.sympy / approximate_degree(D.sympy, self.var, degree)
         return self.__class__(expr, **self.assumptions)
 
-    def approximate_order(self, order):
-        """Reduce order of numerator and denominator to `order`.
+    def approximate_degree(self, degree):
+        """Reduce degree of numerator and denominator to `degree`.
         Warning, this assumes `self.var` is much smaller than 1."""
 
         N, D = self.as_N_D()
 
-        expr = approximate_order(N.sympy, self.var, order) / \
-            approximate_order(D.sympy, self.var, order)
+        expr = approximate_degree(N.sympy, self.var, degree) / \
+            approximate_degree(D.sympy, self.var, degree)
         return self.__class__(expr, **self.assumptions)
 
-    def approximate_taylor(self, var0, degree):
+    def approximate_pade(self, ndegree=2, ddegree=2):
+        """Approximate expression using a Pade series with numerator degree
+        `ndegree` and denominator degree `ddegree`.
+
+        This determines the Pade coeficients from a Taylor series
+        expansion and thus fails if the expression has singularities
+        at the origin.
+
+        """
+
+        expr = approximate_pade(self.sympy, var=self.var,
+                                ndegree=ndegree, ddegree=ddegree)
+        return self.__class__(expr, **self.assumptions)
+
+    def approximate_taylor(self, degree=2, var0=0):
         """Approximate expression using a Taylor series
         around `self.var = var0` to degree `degree`."""
 
         from .expr import expr
 
         var0 = expr(var0)
-        expr = approximate_taylor(self.sympy, self.var, var0.sympy, degree)
+        expr = approximate_taylor(self.sympy, var=self.var,
+                                  degree=degree, var0=var0.sympy)
         return self.__class__(expr, **self.assumptions)
 
-    def approximate(self, method='pade', order=1, numer_order=None, defs=None,
+    def approximate(self, ndegree=2, ddegree=2, method=None, defs=None,
                     threshold=0.01):
         """Approximate an expression.  The order of attempted approximations is:
         1. Fractional powers
         2. Exponential functions
         3. Hyperbolic trigonometric functions
         4. Dominant terms (if `defs` is defined).
 
-        See also `approximate_numer_order`, `approximate_denom_order`, and
-        `approximate_order` methods.
+        See also `approximate_numer_degree`, `approximate_denom_degree`, and
+        `approximate_degree` methods.
         """
 
         terms = self.expand().as_ordered_terms()
         if len(terms) > 1:
             result = 0
             for term in terms:
                 result += term.approximate(method,
-                                           order, numer_order, defs, threshold)
+                                           ndegree, ddegree, defs, threshold)
             return result
 
-        result = self.approximate_fractional_power(method=method,
-                                                   order=order)
-        result = result.approximate_exp(method=method,
-                                        order=order,
-                                        numer_order=numer_order)
-        result = result.approximate_hyperbolic_trig(method=method,
-                                                    order=order,
-                                                    numer_order=numer_order)
+        if method == 'taylor':
+            result = self.approximate_taylor(ndegree)
+        elif method == 'pade':
+            result = self.approximate_pade(ndegree, ddegree)
+        else:
+
+            result = self.approximate_fractional_power(ndegree=ndegree,
+                                                       ddegree=ddegree)
+            result = result.approximate_exp(ndegree=ndegree,
+                                            ddegree=ddegree)
+            result = result.approximate_hyperbolic_trig(ndegree=ndegree,
+                                                        ddegree=ddegree)
         if defs is not None:
             result = result.approximate_dominant(defs, threshold)
         return result