"Fossies" - the Fresh Open Source Software Archive

Member "dmelt/python/packages/sympy/mpmath/calculus/approximation.py" (10 Apr 2019, 8743 Bytes) of package /linux/misc/dmelt-2.4.zip:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "approximation.py" see the Fossies "Dox" file reference documentation.

    1 from calculus import defun
    2 
    3 #----------------------------------------------------------------------------#
    4 #                              Approximation methods                         #
    5 #----------------------------------------------------------------------------#
    6 
    7 # The Chebyshev approximation formula is given at:
    8 # http://mathworld.wolfram.com/ChebyshevApproximationFormula.html
    9 
   10 # The only major changes in the following code is that we return the
   11 # expanded polynomial coefficients instead of Chebyshev coefficients,
   12 # and that we automatically transform [a,b] -> [-1,1] and back
   13 # for convenience.
   14 
   15 # Coefficient in Chebyshev approximation
   16 def chebcoeff(ctx,f,a,b,j,N):
   17     s = ctx.mpf(0)
   18     h = ctx.mpf(0.5)
   19     for k in range(1, N+1):
   20         t = ctx.cos(ctx.pi*(k-h)/N)
   21         s += f(t*(b-a)*h + (b+a)*h) * ctx.cos(ctx.pi*j*(k-h)/N)
   22     return 2*s/N
   23 
   24 # Generate Chebyshev polynomials T_n(ax+b) in expanded form
   25 def chebT(ctx, a=1, b=0):
   26     Tb = [1]
   27     yield Tb
   28     Ta = [b, a]
   29     while 1:
   30         yield Ta
   31         # Recurrence: T[n+1](ax+b) = 2*(ax+b)*T[n](ax+b) - T[n-1](ax+b)
   32         Tmp = [0] + [2*a*t for t in Ta]
   33         for i, c in enumerate(Ta): Tmp[i] += 2*b*c
   34         for i, c in enumerate(Tb): Tmp[i] -= c
   35         Ta, Tb = Tmp, Ta
   36 
   37 @defun
   38 def chebyfit(ctx, f, interval, N, error=False):
   39     r"""
   40     Computes a polynomial of degree `N-1` that approximates the
   41     given function `f` on the interval `[a, b]`. With ``error=True``,
   42     :func:`chebyfit` also returns an accurate estimate of the
   43     maximum absolute error; that is, the maximum value of
   44     `|f(x) - P(x)|` for `x \in [a, b]`.
   45 
   46     :func:`chebyfit` uses the Chebyshev approximation formula,
   47     which gives a nearly optimal solution: that is, the maximum
   48     error of the approximating polynomial is very close to
   49     the smallest possible for any polynomial of the same degree.
   50 
   51     Chebyshev approximation is very useful if one needs repeated
   52     evaluation of an expensive function, such as function defined
   53     implicitly by an integral or a differential equation. (For
   54     example, it could be used to turn a slow mpmath function
   55     into a fast machine-precision version of the same.)
   56 
   57     **Examples**
   58 
   59     Here we use :func:`chebyfit` to generate a low-degree approximation
   60     of `f(x) = \cos(x)`, valid on the interval `[1, 2]`::
   61 
   62         >>> from mpmath import *
   63         >>> mp.dps = 15; mp.pretty = True
   64         >>> poly, err = chebyfit(cos, [1, 2], 5, error=True)
   65         >>> nprint(poly)
   66         [0.00291682, 0.146166, -0.732491, 0.174141, 0.949553]
   67         >>> nprint(err, 12)
   68         1.61351758081e-5
   69 
   70     The polynomial can be evaluated using ``polyval``::
   71 
   72         >>> nprint(polyval(poly, 1.6), 12)
   73         -0.0291858904138
   74         >>> nprint(cos(1.6), 12)
   75         -0.0291995223013
   76 
   77     Sampling the true error at 1000 points shows that the error
   78     estimate generated by ``chebyfit`` is remarkably good::
   79 
   80         >>> error = lambda x: abs(cos(x) - polyval(poly, x))
   81         >>> nprint(max([error(1+n/1000.) for n in range(1000)]), 12)
   82         1.61349954245e-5
   83 
   84     **Choice of degree**
   85 
   86     The degree `N` can be set arbitrarily high, to obtain an
   87     arbitrarily good approximation. As a rule of thumb, an
   88     `N`-term Chebyshev approximation is good to `N/(b-a)` decimal
   89     places on a unit interval (although this depends on how
   90     well-behaved `f` is). The cost grows accordingly: ``chebyfit``
   91     evaluates the function `(N^2)/2` times to compute the
   92     coefficients and an additional `N` times to estimate the error.
   93 
   94     **Possible issues**
   95 
   96     One should be careful to use a sufficiently high working
   97     precision both when calling ``chebyfit`` and when evaluating
   98     the resulting polynomial, as the polynomial is sometimes
   99     ill-conditioned. It is for example difficult to reach
  100     15-digit accuracy when evaluating the polynomial using
  101     machine precision floats, no matter the theoretical
  102     accuracy of the polynomial. (The option to return the
  103     coefficients in Chebyshev form should be made available
  104     in the future.)
  105 
  106     It is important to note the Chebyshev approximation works
  107     poorly if `f` is not smooth. A function containing singularities,
  108     rapid oscillation, etc can be approximated more effectively by
  109     multiplying it by a weight function that cancels out the
  110     nonsmooth features, or by dividing the interval into several
  111     segments.
  112     """
  113     a, b = ctx._as_points(interval)
  114     orig = ctx.prec
  115     try:
  116         ctx.prec = orig + int(N**0.5) + 20
  117         c = [chebcoeff(ctx,f,a,b,k,N) for k in range(N)]
  118         d = [ctx.zero] * N
  119         d[0] = -c[0]/2
  120         h = ctx.mpf(0.5)
  121         T = chebT(ctx, ctx.mpf(2)/(b-a), ctx.mpf(-1)*(b+a)/(b-a))
  122         for k in range(N):
  123             Tk = T.next()
  124             for i in range(len(Tk)):
  125                 d[i] += c[k]*Tk[i]
  126         d = d[::-1]
  127         # Estimate maximum error
  128         err = ctx.zero
  129         for k in range(N):
  130             x = ctx.cos(ctx.pi*k/N) * (b-a)*h + (b+a)*h
  131             err = max(err, abs(f(x) - ctx.polyval(d, x)))
  132     finally:
  133         ctx.prec = orig
  134     if error:
  135         return d, +err
  136     else:
  137         return d
  138 
  139 @defun
  140 def fourier(ctx, f, interval, N):
  141     r"""
  142     Computes the Fourier series of degree `N` of the given function
  143     on the interval `[a, b]`. More precisely, :func:`fourier` returns
  144     two lists `(c, s)` of coefficients (the cosine series and sine
  145     series, respectively), such that
  146 
  147     .. math ::
  148 
  149         f(x) \sim \sum_{k=0}^N
  150             c_k \cos(k m) + s_k \sin(k m)
  151 
  152     where `m = 2 \pi / (b-a)`.
  153 
  154     Note that many texts define the first coefficient as `2 c_0` instead
  155     of `c_0`. The easiest way to evaluate the computed series correctly
  156     is to pass it to :func:`fourierval`.
  157 
  158     **Examples**
  159 
  160     The function `f(x) = x` has a simple Fourier series on the standard
  161     interval `[-\pi, \pi]`. The cosine coefficients are all zero (because
  162     the function has odd symmetry), and the sine coefficients are
  163     rational numbers::
  164 
  165         >>> from mpmath import *
  166         >>> mp.dps = 15; mp.pretty = True
  167         >>> c, s = fourier(lambda x: x, [-pi, pi], 5)
  168         >>> nprint(c)
  169         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  170         >>> nprint(s)
  171         [0.0, 2.0, -1.0, 0.666667, -0.5, 0.4]
  172 
  173     This computes a Fourier series of a nonsymmetric function on
  174     a nonstandard interval::
  175 
  176         >>> I = [-1, 1.5]
  177         >>> f = lambda x: x**2 - 4*x + 1
  178         >>> cs = fourier(f, I, 4)
  179         >>> nprint(cs[0])
  180         [0.583333, 1.12479, -1.27552, 0.904708, -0.441296]
  181         >>> nprint(cs[1])
  182         [0.0, -2.6255, 0.580905, 0.219974, -0.540057]
  183 
  184     It is instructive to plot a function along with its truncated
  185     Fourier series::
  186 
  187         >>> plot([f, lambda x: fourierval(cs, I, x)], I) #doctest: +SKIP
  188 
  189     Fourier series generally converge slowly (and may not converge
  190     pointwise). For example, if `f(x) = \cosh(x)`, a 10-term Fourier
  191     series gives an `L^2` error corresponding to 2-digit accuracy::
  192 
  193         >>> I = [-1, 1]
  194         >>> cs = fourier(cosh, I, 9)
  195         >>> g = lambda x: (cosh(x) - fourierval(cs, I, x))**2
  196         >>> nprint(sqrt(quad(g, I)))
  197         0.00467963
  198 
  199     :func:`fourier` uses numerical quadrature. For nonsmooth functions,
  200     the accuracy (and speed) can be improved by including all singular
  201     points in the interval specification::
  202 
  203         >>> nprint(fourier(abs, [-1, 1], 0), 10)
  204         ([0.5000441648], [0.0])
  205         >>> nprint(fourier(abs, [-1, 0, 1], 0), 10)
  206         ([0.5], [0.0])
  207 
  208     """
  209     interval = ctx._as_points(interval)
  210     a = interval[0]
  211     b = interval[-1]
  212     L = b-a
  213     cos_series = []
  214     sin_series = []
  215     cutoff = ctx.eps*10
  216     for n in xrange(N+1):
  217         m = 2*n*ctx.pi/L
  218         an = 2*ctx.quadgl(lambda t: f(t)*ctx.cos(m*t), interval)/L
  219         bn = 2*ctx.quadgl(lambda t: f(t)*ctx.sin(m*t), interval)/L
  220         if n == 0:
  221             an /= 2
  222         if abs(an) < cutoff: an = ctx.zero
  223         if abs(bn) < cutoff: bn = ctx.zero
  224         cos_series.append(an)
  225         sin_series.append(bn)
  226     return cos_series, sin_series
  227 
  228 @defun
  229 def fourierval(ctx, series, interval, x):
  230     """
  231     Evaluates a Fourier series (in the format computed by
  232     by :func:`fourier` for the given interval) at the point `x`.
  233 
  234     The series should be a pair `(c, s)` where `c` is the
  235     cosine series and `s` is the sine series. The two lists
  236     need not have the same length.
  237     """
  238     cs, ss = series
  239     ab = ctx._as_points(interval)
  240     a = interval[0]
  241     b = interval[-1]
  242     m = 2*ctx.pi/(ab[-1]-ab[0])
  243     s = ctx.zero
  244     s += ctx.fsum(cs[n]*ctx.cos(m*n*x) for n in xrange(len(cs)) if cs[n])
  245     s += ctx.fsum(ss[n]*ctx.sin(m*n*x) for n in xrange(len(ss)) if ss[n])
  246     return s