Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/approximation.py: 14%
73 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from ..libmp.backend import xrange
2from .calculus import defun
4#----------------------------------------------------------------------------#
5# Approximation methods #
6#----------------------------------------------------------------------------#
8# The Chebyshev approximation formula is given at:
9# http://mathworld.wolfram.com/ChebyshevApproximationFormula.html
11# The only major changes in the following code is that we return the
12# expanded polynomial coefficients instead of Chebyshev coefficients,
13# and that we automatically transform [a,b] -> [-1,1] and back
14# for convenience.
16# Coefficient in Chebyshev approximation
17def chebcoeff(ctx,f,a,b,j,N):
18 s = ctx.mpf(0)
19 h = ctx.mpf(0.5)
20 for k in range(1, N+1):
21 t = ctx.cospi((k-h)/N)
22 s += f(t*(b-a)*h + (b+a)*h) * ctx.cospi(j*(k-h)/N)
23 return 2*s/N
25# Generate Chebyshev polynomials T_n(ax+b) in expanded form
26def chebT(ctx, a=1, b=0):
27 Tb = [1]
28 yield Tb
29 Ta = [b, a]
30 while 1:
31 yield Ta
32 # Recurrence: T[n+1](ax+b) = 2*(ax+b)*T[n](ax+b) - T[n-1](ax+b)
33 Tmp = [0] + [2*a*t for t in Ta]
34 for i, c in enumerate(Ta): Tmp[i] += 2*b*c
35 for i, c in enumerate(Tb): Tmp[i] -= c
36 Ta, Tb = Tmp, Ta
38@defun
39def chebyfit(ctx, f, interval, N, error=False):
40 r"""
41 Computes a polynomial of degree `N-1` that approximates the
42 given function `f` on the interval `[a, b]`. With ``error=True``,
43 :func:`~mpmath.chebyfit` also returns an accurate estimate of the
44 maximum absolute error; that is, the maximum value of
45 `|f(x) - P(x)|` for `x \in [a, b]`.
47 :func:`~mpmath.chebyfit` uses the Chebyshev approximation formula,
48 which gives a nearly optimal solution: that is, the maximum
49 error of the approximating polynomial is very close to
50 the smallest possible for any polynomial of the same degree.
52 Chebyshev approximation is very useful if one needs repeated
53 evaluation of an expensive function, such as function defined
54 implicitly by an integral or a differential equation. (For
55 example, it could be used to turn a slow mpmath function
56 into a fast machine-precision version of the same.)
58 **Examples**
60 Here we use :func:`~mpmath.chebyfit` to generate a low-degree approximation
61 of `f(x) = \cos(x)`, valid on the interval `[1, 2]`::
63 >>> from mpmath import *
64 >>> mp.dps = 15; mp.pretty = True
65 >>> poly, err = chebyfit(cos, [1, 2], 5, error=True)
66 >>> nprint(poly)
67 [0.00291682, 0.146166, -0.732491, 0.174141, 0.949553]
68 >>> nprint(err, 12)
69 1.61351758081e-5
71 The polynomial can be evaluated using ``polyval``::
73 >>> nprint(polyval(poly, 1.6), 12)
74 -0.0291858904138
75 >>> nprint(cos(1.6), 12)
76 -0.0291995223013
78 Sampling the true error at 1000 points shows that the error
79 estimate generated by ``chebyfit`` is remarkably good::
81 >>> error = lambda x: abs(cos(x) - polyval(poly, x))
82 >>> nprint(max([error(1+n/1000.) for n in range(1000)]), 12)
83 1.61349954245e-5
85 **Choice of degree**
87 The degree `N` can be set arbitrarily high, to obtain an
88 arbitrarily good approximation. As a rule of thumb, an
89 `N`-term Chebyshev approximation is good to `N/(b-a)` decimal
90 places on a unit interval (although this depends on how
91 well-behaved `f` is). The cost grows accordingly: ``chebyfit``
92 evaluates the function `(N^2)/2` times to compute the
93 coefficients and an additional `N` times to estimate the error.
95 **Possible issues**
97 One should be careful to use a sufficiently high working
98 precision both when calling ``chebyfit`` and when evaluating
99 the resulting polynomial, as the polynomial is sometimes
100 ill-conditioned. It is for example difficult to reach
101 15-digit accuracy when evaluating the polynomial using
102 machine precision floats, no matter the theoretical
103 accuracy of the polynomial. (The option to return the
104 coefficients in Chebyshev form should be made available
105 in the future.)
107 It is important to note the Chebyshev approximation works
108 poorly if `f` is not smooth. A function containing singularities,
109 rapid oscillation, etc can be approximated more effectively by
110 multiplying it by a weight function that cancels out the
111 nonsmooth features, or by dividing the interval into several
112 segments.
113 """
114 a, b = ctx._as_points(interval)
115 orig = ctx.prec
116 try:
117 ctx.prec = orig + int(N**0.5) + 20
118 c = [chebcoeff(ctx,f,a,b,k,N) for k in range(N)]
119 d = [ctx.zero] * N
120 d[0] = -c[0]/2
121 h = ctx.mpf(0.5)
122 T = chebT(ctx, ctx.mpf(2)/(b-a), ctx.mpf(-1)*(b+a)/(b-a))
123 for (k, Tk) in zip(range(N), T):
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
139@defun
140def 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:`~mpmath.fourier` returns
144 two lists `(c, s)` of coefficients (the cosine series and sine
145 series, respectively), such that
147 .. math ::
149 f(x) \sim \sum_{k=0}^N
150 c_k \cos(k m x) + s_k \sin(k m x)
152 where `m = 2 \pi / (b-a)`.
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:`~mpmath.fourierval`.
158 **Examples**
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::
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]
173 This computes a Fourier series of a nonsymmetric function on
174 a nonstandard interval::
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]
184 It is instructive to plot a function along with its truncated
185 Fourier series::
187 >>> plot([f, lambda x: fourierval(cs, I, x)], I) #doctest: +SKIP
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::
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
199 :func:`~mpmath.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::
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])
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
228@defun
229def fourierval(ctx, series, interval, x):
230 """
231 Evaluates a Fourier series (in the format computed by
232 by :func:`~mpmath.fourier` for the given interval) at the point `x`.
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