## "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
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