# Lesson 33. Approximation of functions
# --------------------------------------
#
# Suppose we have to work with a certain function f(x), defined on an interval. It is a well-behaved "smooth"
# function, but is not easy to compute from its definition. For example, it might require an integration or solving
# an implicit equation. Since we'll need to compute its values many times, we are willing to put some effort into
# developing an efficient way of obtaining numerical approximations for f(x).
#
# We'll consider the following function on the interval [-1, 1]. Of course Maple can evaluate it quite efficiently,
# but we'll pretend that it can't: for a more realistic example, some of our computations would take a
# considerable amount of computer time.
> restart: f:= x -> arctan(2*x+1):
> plot(f(x), x=-1..1);
** Maple V Graphics **
# The goal is to find an approximation of f(x) by a polynomial with an error of at most 10^(-8), requiring as few
# arithmetic operations as possible.
#
# Hearing the words "polynomial approximation", the Calculus student will first think of Taylor series. However,
# these turn out not to be useful for our purpose.
> taylor(f(x), x=0, 11);
2 3 5 6 7 9 10
1/4 Pi + x - x + 2/3 x - 4/5 x + 4/3 x - 8/7 x + 16/9 x - 16/5 x
11
+ O(x )
> TaylorApp:= unapply(convert(",polynom), x);
TaylorApp :=
2 3 5 6 7 9 10
x -> 1/4 Pi + x - x + 2/3 x - 4/5 x + 4/3 x - 8/7 x + 16/9 x - 16/5 x
> plot(f(x)-TaylorApp(x), x = -1 .. 1);
** Maple V Graphics **
# The approximation is very good for x near 0, but poor near the ends of the interval, so the maximum error is
# about 2.8. In fact, the radius of convergence of the Taylor series is sqrt(2)/2, and for |x| larger than this the
# Taylor series is basically useless.
> MaxTaylorError:= evalf(f(-1)-TaylorApp(-1));
MaxTaylorError := 2.797457641
# The next type of series we looked at was a Fourier series.
# For convenience we use a change of variables to get a function defined on [0,Pi].
> fp:= t -> f(-1 + 2/Pi * t):
# Recall that for a Fourier series we consider our function fp(t) to be made into a periodic function. The best
# way to do that here is to make it an even function. Otherwise it would not be continuous, and the Fourier series
# would not converge to the value of the functon at the discontinuity.
> plots[display]({plot(fp(abs(t)), t = -Pi .. Pi),\
plot(fp(t+2*Pi), t=-2*Pi .. -Pi),\
plot(fp(-t+2*Pi), t=Pi .. 2*Pi)});
** Maple V Graphics **
# Since it's an even function, the sin terms are 0 by symmetry, and we only need to integrate over [0,Pi] to get the
# cos terms.
> for kk from 0 to 10 do \
a[kk]:= evalf(2/Pi * Int(fp(t)*cos(kk*t), t=0..Pi)) od:\
FourierApp:= unapply(a[0]/2 + sum(a[k]*cos(k*(x+1)*Pi/2), k = 1 .. 10), x);
FourierApp := x -> .5392550495 - .8754281782 cos(1/2 (x + 1) Pi)
- .2613646026 cos((x + 1) Pi) - .09024319630 cos(3/2 (x + 1) Pi)
- .01778157084 cos(2 (x + 1) Pi) - .01262891368 cos(5/2 (x + 1) Pi)
- .005502151118 cos(3 (x + 1) Pi) - .008769972060 cos(7/2 (x + 1) Pi)
- .004907326798 cos(4 (x + 1) Pi) - .006019734224 cos(9/2 (x + 1) Pi)
- .003256263088 cos(5 (x + 1) Pi)
> plot(f(x) - FourierApp(x), x=-1 .. 1);
** Maple V Graphics **
# This time the approximation is quite good. The maximum error is again at x=-1.
> MaxFourierError:= evalf(FourierApp(-1)-f(-1));
MaxFourierError := .0387513041
# Recall from Lessons 25 and 26 that the error in the n'th partial sum of the Fourier series is typically O(1/n) at a
# sharp "corner" where the derivative of the function is discontinuous, but converges faster when the function is
# smooth. A different change of variables could "smooth out" the graph at x=-1 and x=1. The one we want is x
# = cos(t), 0 <= t <= Pi.
# In fact, f(cos(t)) is an analytic function of t. There is an additional, very important, benefit: cos(k t) is a
# polynomial in cos(t), so our approximation becomes a polynomial in x rather than a trigonometric function.
> seq(expand(cos(kk*t)), kk = 2 .. 10);
2 3 4 2
2 cos(t) - 1, 4 cos(t) - 3 cos(t), 8 cos(t) - 8 cos(t) + 1,
5 3
16 cos(t) - 20 cos(t) + 5 cos(t),
6 4 2
32 cos(t) - 48 cos(t) + 18 cos(t) - 1,
7 5 3
64 cos(t) - 112 cos(t) + 56 cos(t) - 7 cos(t),
8 6 4 2
128 cos(t) - 256 cos(t) + 160 cos(t) - 32 cos(t) + 1,
9 7 5 3
256 cos(t) - 576 cos(t) + 432 cos(t) - 120 cos(t) + 9 cos(t),
10 8 6 4 2
512 cos(t) - 1280 cos(t) + 1120 cos(t) - 400 cos(t) + 50 cos(t) - 1
# The functions expressing cos(k t) as a polynomial in cos(t) are called "Chebyshev polynomials of the first kind",
# and the "Fourier" series of f(x) as a sum of Chebyshev polynomials is called a Chebyshev series.
> fp:= unapply(f(cos(t)), t);
fp := t -> arctan(2 cos(t) + 1)
> plot(fp(t), t = -2*Pi .. 2*Pi);
** Maple V Graphics **
> for kk from 0 to 10 do \
a[kk]:= evalf(2/Pi * Int(fp(t)*cos(kk*t), t=0..Pi)) od:\
a[0]/2 + sum(a[k]*cos(k*t), k = 1 .. 10);
.4522784470 + 1.058171027 cos(t) - .2720196496 cos(2 t) - .02881149458 cos(3 t)
+ .05817102726 cos(4 t) - .01794454247 cos(5 t) - .005730457514 cos(6 t)
+ .006960372170 cos(7 t) - .001644114297 cos(8 t) - .001121821048 cos(9 t)
+ .0009741971124 cos(10 t)
> ChebApp:= unapply(subs(cos(t)=x, expand(")),x);
3 4 2
ChebApp := x -> .7680442384 x + .9960638042 x + .0876930463 x - 1.011234240 x
7 8 6
+ 1.091632743 x - 1.457418934 x + 1.328619386 x + .7855812700
10 9 5
+ .4987889215 x - .2871861883 x - 1.551301056 x
> plot(f(x)-ChebApp(x), x=-1 .. 1);
** Maple V Graphics **
# The maximum error here is much smaller: about .0004. Moreover the errors are spread out over the interval,
# rather than being concentrated in one part of it.
#
# Maple has its own built-in facilities for Chebyshev (and other) approximations to functions, in the
# "numapprox" package.
> with(numapprox):
# We'll be trying to get some very accurate approximations to f(x) now, so to reduce roundoff error I'll increase
# Digits.
> Digits:= 15:
# The "numapprox" commands require a function that always returns numerical values, so I'll make a version of
# f that does that.
> ff:= x -> evalf(f(x)):
# The "chebyshev" command produces an approximation of a function over a given interval in terms of
# Chebyshev polynomials. The first input is the function, the second specifies the interval, and the third is the
# desired accuracy (call it "eps"). The result is a sum of Chebyshev polynomials (written as T(k, x)) times
# floating-point coefficients, with only the terms in the series having coefficients of absolute value at least eps/5
# times the largest coefficient.
> chebs:=chebyshev(ff,x = -1 .. 1, 1E-7);
chebs := .452278446998695 T(0, x) + 1.05817102767564 T(1, x)
- .272019649701265 T(2, x) - .0288114956717208 T(3, x)
+ .0581710297420256 T(4, x) - .0179445429927063 T(5, x)
- .00573046437003628 T(6, x) + .00696038436235612 T(7, x)
- .00164411184515822 T(8, x) - .00112186268711132 T(9, x)
+ .000974253735967679 T(10, x) - .000137699438193821 T(11, x)
- .000216232125443209 T(12, x) + .000143076436570987 T(13, x)
-5
- .374270813849617*10 T(14, x) - .0000422249399142610 T(15, x)
-5
+ .0000218511788069392 T(16, x) + .530515292994653*10 T(17, x)
- .0000151528101015989 T(18, x)
# To see these Chebyshev polynomials as polynomials, you substitute "orthopoly[T]" for "T", i.e. the Chebyshev
# polynomials themselves are defined as T(k,x) in the "orthopoly" package.
> seq(orthopoly[T](kk,x) = subs(cos(t)=x, expand(cos(kk*t))), kk = 0 .. 5);
2 2 3 3
1 = 1, x = x, 2 x - 1 = 2 x - 1, 4 x - 3 x = 4 x - 3 x,
4 2 4 2 5 3 5 3
8 x - 8 x + 1 = 8 x - 8 x + 1, 16 x - 20 x + 5 x = 16 x - 20 x + 5 x
> ChebApp2:= unapply(eval(subs(T=orthopoly[T],chebs)),x);
18 15
ChebApp2 := x -> - 1.98610912563677 x - 2.16944705082740 x
14 13 12
- 19.6525337212140 x + 5.76620025425322 x + 21.2625872893756 x
11 17 16
- 8.33824286815898 x + .347678502416976 x + 9.65351049251125 x
3 4
+ .785405739802500 + .657701253455242 x + .0361008273463846 x
2 7 8
- 1.00128140864535 x - 2.32160994483594 x + 2.45079940402496 x
6 10 9
+ .926501313534097 x - 12.2431645830034 x + 6.72383384987049 x
5
- .649053622476610 x + 1.00016159420086 x
> plot(f(x)-ChebApp2(x),x=-1..1);
** Maple V Graphics **
# The maximum error here is approximately 1.2 * 10^(-5). Rather disappointing: I would have hoped for
# something more like 10^(-7).
> chebs2:= chebyshev(ff,x=-1..1,1E-8);
chebs2 := .452278447151191 T(0, x) + 1.05817102727149 T(1, x)
- .272019649514068 T(2, x) - .0288114945860670 T(3, x)
+ .0581710272714922 T(4, x) - .0179445424708701 T(5, x)
- .00573045751510577 T(6, x) + .00696037216927649 T(7, x)
- .00164411429807496 T(8, x) - .00112182104832470 T(9, x)
+ .000974197112900721 T(10, x) - .000137741903096464 T(11, x)
- .000215989116121855 T(12, x) + .000142833053852068 T(13, x)
-5
- .413079194212294*10 T(14, x) - .0000408504603295013 T(15, x)
-5
+ .0000209405094970248 T(16, x) + .233904648046266*10 T(17, x)
-5 -5
- .757640505947170*10 T(18, x) + .296610644775340*10 T(19, x)
-6 -5
+ .910669348962297*10 T(20, x) - .137447961974368*10 T(21, x)
-6 -6
+ .388083734916542*10 T(22, x) + .243382924191194*10 T(23, x)
-6 -7
- .243009422309937*10 T(24, x) + .424644369389531*10 T(25, x)
-7 -7
+ .566240774495625*10 T(26, x) - .416389150745983*10 T(27, x)
-8 -7
+ .245005884385748*10 T(28, x) + .121978131685512*10 T(29, x)
-8 -8
- .685347846457773*10 T(30, x) + .249134437725203*10 T(32, x)
# It seems that "chebyshev" didn't really perform as advertised: I don't know why the first time it left out the
# terms in T(19,x) to T(27,x) which all have coefficients greater than eps/5 times the largest coefficient.
# Somehow it's deciding where to stop looking at more coefficients, using some sort of estimate of the error, and
# the estimate was over-optimistic in this case. Let's take the polynomial approximation corresponding to
# "chebs2". One added wrinkle: it's more efficient (i.e. requires fewer arithmetic operations) to evaluate a
# polynomial in "Horner form". The "hornerform" command in "numapprox" converts a polynomial to that form.
> hornerform(2+3*x+4*x^2+5*x^3);
2 + (3 + (4 + 5 x) x) x
> ChebApp3:= unapply(hornerform(eval(subs(T=orthopoly[T],chebs2))),x);
ChebApp3 := x -> .785398163378253 + (1.00000003708117 + ( - 1.00000001538419 + (
.666659195615746 + (.45105973175186*10^(-5) + ( - .799557318661753 + (
1.33302399177422 + ( - 1.15502600715574 + (.00948684558167 + (
1.96627156266228 + ( - 3.36275456728765 + (1.08405204349057 + (
1.74944162903890 + (6.78802630133392 + ( - 3.41995388010094 + (
- 59.7933980014170 + (62.4078750110457 + (168.294379003530 + (
- 245.384183110292 + ( - 274.152429280531 + (519.592552650623 + (
292.036755683805 + ( - 712.517408172944 + ( - 208.385163855401 + (
668.685711468569 + (96.7255275192457 + ( - 430.595196927034 + (
- 26.5332004545944 + (183.078108622744
2
+ (3.27432554010285 + ( - 46.4804037271340 + 5.35012131168548 x ) x) x) x)
x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x) x
) x
# This one is going to be so accurate that I'll multiply the error by 10^8 before plotting (otherwise Maple
# wouldn't label the scale on the y axis).
> plot((f(x)-ChebApp3(x))*1E8, x=-1..1);
** Maple V Graphics **
# The maximum error is about 1.7*10^(-9).
#
# Sometimes it's better to use an approximation by a rational function rather than by a polynomial. It's actually
# not much help in the case of our f, so let's try a different example, still on the interval [-1,1].
> g:= x -> 1/(exp(3*x) + 1 + x):
> plot(g(x), x=-1 .. 1);
** Maple V Graphics **
# The left endpoint is very close to a vertical asymptote of the function. We don't expect a Taylor series to do
# very well near there.
> convert(taylor(g(x), x=0, 10),polynom);
2 3 11 4 121 5 419 6 527 7 5473 8
1/2 - x + 7/8 x - 5/8 x + ---- x - --- x + --- x - --- x + ---- x
16 160 640 896 8960
191 9
- --- x
320
> TayApp:= unapply(evalf("),x);
2
TayApp := x -> .500000000000000 - 1. x + .875000000000000 x
3 4 5
- .625000000000000 x + .687500000000000 x - .756250000000000 x
6 7 8
+ .654687500000000 x - .588169642857143 x + .610825892857143 x
9
- .596875000000000 x
# The (m,n) Pade approximation of f at x=0 is a rational function p(x)/q(x), with p a polynomial of degree <= m
# and q a polynomial of degree <= n, whose Taylor series at x=0 agrees as much as possible with the Taylor series
# of f. In normal cases, the series agree through the term of degree m+n. The "pade" command in "numapprox"
# will calculate it.
> p:=pade(g(x),x=0,[5,4]);
28 2 15 3 27 4 27 5
1/2 - ---- x + 3/8 x - --- x + ---- x + ----- x
41 164 9184 11480
p := ---------------------------------------------------
26 11 2 81 3 267 4
1 + ---- x + ---- x + --- x - ---- x
41 41 164 4592
> taylor(",x=0,11);
2 3 11 4 121 5 419 6 527 7 5473 8
1/2 - x + 7/8 x - 5/8 x + ---- x - --- x + --- x - --- x + ---- x
16 160 640 896 8960
191 9 1596373 10 11
- --- x + ------- x + O(x )
320 2938880
> PadeApp:= unapply(evalf(p),x);
2
PadeApp := x -> (.500000000000000 - .682926829268293 x + .375000000000000 x
3 4 5
- .0914634146341463 x + .00293989547038328 x + .00235191637630662 x )
/ 2 3
/ (1. + .634146341463415 x + .268292682926829 x + .493902439024390 x
/
4
- .0581445993031359 x )
> plot(g(x)-TayApp(x), x=-1.. 1);
** Maple V Graphics **
> plot(g(x)-PadeApp(x), x=-1 .. 1);
** Maple V Graphics **
# Much better!
> gg:= x -> evalf(g(x)):
# Now to try Chebyshev. Since the function is so badly behaved, I won't be too demanding in terms of accuracy.
> ch:= chebyshev(gg, x=-1 .. 1, 1E-4);
ch := 2.65190894460027 T(0, 1.00000000000000 x)
- 4.53750217235296 T(1, 1.00000000000000 x)
+ 3.32124279412398 T(2, 1.00000000000000 x)
- 2.43058149985833 T(3, 1.00000000000000 x)
+ 1.81989826687969 T(4, 1.00000000000000 x)
- 1.35963393371419 T(5, 1.00000000000000 x)
+ 1.01149398103648 T(6, 1.00000000000000 x)
- .753146988669057 T(7, 1.00000000000000 x)
+ .561205994081705 T(8, 1.00000000000000 x)
- .418077043436990 T(9, 1.00000000000000 x)
+ .311412815862561 T(10, 1.00000000000000 x)
- .231976434906242 T(11, 1.00000000000000 x)
+ .172806151325234 T(12, 1.00000000000000 x)
- .128726591050350 T(13, 1.00000000000000 x)
+ .0958906804211621 T(14, 1.00000000000000 x)
- .0714308599824060 T(15, 1.00000000000000 x)
+ .0532102630020011 T(16, 1.00000000000000 x)
- .0396373542776992 T(17, 1.00000000000000 x)
+ .0295266332707976 T(18, 1.00000000000000 x)
- .0219949639498650 T(19, 1.00000000000000 x)
+ .0163844765057379 T(20, 1.00000000000000 x)
- .0122051150093616 T(21, 1.00000000000000 x)
+ .00909182746859014 T(22, 1.00000000000000 x)
- .00677267909042103 T(23, 1.00000000000000 x)
+ .00504510037566792 T(24, 1.00000000000000 x)
- .00375819341783964 T(25, 1.00000000000000 x)
+ .00279955142861861 T(26, 1.00000000000000 x)
- .00208544040272035 T(27, 1.00000000000000 x)
+ .00155348525524814 T(28, 1.00000000000000 x)
- .00115722157672665 T(29, 1.00000000000000 x)
+ .000862037129665312 T(30, 1.00000000000000 x)
- .000642148591602401 T(31, 1.00000000000000 x)
+ .000478349472721025 T(32, 1.00000000000000 x)
- .000356332496544850 T(33, 1.00000000000000 x)
+ .000265439909691009 T(34, 1.00000000000000 x)
- .000197732564629368 T(35, 1.00000000000000 x)
+ .000147296463091237 T(36, 1.00000000000000 x)
- .000109726189136820 T(37, 1.00000000000000 x)
> ChebApp:= unapply(eval(subs(T=orthopoly[T],ch)),x);
2 3 5
ChebApp := x -> .8384829447599 x - .259658788547586 x - 27.248645487571 x
6 7 9
- 190.752459951669 x + 901.65588122390 x - 17573.7886392227 x
10 29
- 66739.7006124610 x + .50005256730982 - .110573224245382*10^10 x
31 28
+ .660269114076604*10^9 x + .650674953880608*10^9 x
30 32
- .403333325308618*10^9 x + .169538447303325*10^9 x
20 22
+ .202958370133194*10^9 x - .415517297800262*10^9 x
24 26
+ .644779312942766*10^9 x - .752959065958479*10^9 x
34 36
- .432694999467679*10^8 x + 5061067.93434667 x
19 21
+ .165651731561531*10^9 x - .419328800294001*10^9 x
23 25
+ .810527319754546*10^9 x - .119515475888133*10^10 x
37 33
- 7540326.30171764 x - .268235592710130*10^9 x
35 27
+ .663509986971388*10^8 x + .133335036203000*10^10 x
4 8 14
+ 4.8918784716771 x + 4586.86569789265 x - 4306307.58651251 x
18 15
- .751013304351820*10^8 x + .111803871416771*10^8 x
13 12 11
- 1849490.97494728 x + 641850.66588589 x + 218259.306780854 x
17 16
- .496902703297638*10^8 x + .208951715978786*10^8 x
- 1.00149928091854 x
# Note the large coefficients. Even just calculating this is going to be a problem when x is near 1 or -1.
#
# The (m,n) Chebyshev-Pade approximation of g on an interval [a,b] is a rational function p(x)/q(x) with p a
# polynomial of degree at most m and q a polynomial of degree at most n, expressed as combinations of
# Chebyshev polynomials, such that the Chebyshev series of p(x)/q(x) agrees as far as possible with the
# Chebyshev series of g(x) on the interval [a,b]. The "chebpade" command calculates it.
> chebpade(gg,x=-1 .. 1, [10,10]);
(.585405067833567 T(0, 1.00000000000000 x)
- .626691609109453 T(1, 1.00000000000000 x)
+ .165346665761609 T(2, 1.00000000000000 x)
- .0246412736504998 T(3, 1.00000000000000 x)
+ .00203179497996430 T(4, 1.00000000000000 x)
/
- .0000728073888489652 T(5, 1.00000000000000 x)) / (
/
T(0, 1.00000000000000 x) + .955353160300308 T(1, 1.00000000000000 x)
+ .155133407851009 T(2, 1.00000000000000 x)
+ .127396742723925 T(3, 1.00000000000000 x)
- .000686194849958047 T(4, 1.00000000000000 x)
+ .00178684524866449 T(5, 1.00000000000000 x))
# I'm surprised that this only involves terms of degree 5 in the numerator and denominator, when terms of
# degree up to 10 were allowed. Numerical difficulties may be involved.
> ChebPade:= unapply(eval(subs(T=orthopoly[T],")),x);
2
ChebPade := x -> (.422090197051922 - .553131825102199 x + .314438971683504 x
3 4 5 /
- .0971089468250199 x + .0162543598397144 x - .00116491822158344 x ) /
/
2
(.844180397299033 + .582097158371855 x + .315756374501682 x
3 4 5
+ .473850065922410 x - .00548955879966438 x + .0285895239786318 x )
> plot(gg(x)-ChebPade(x), x=-1 .. 1);
** Maple V Graphics **
> gg(-1)-ChebPade(-1);
-7
.270432*10
>
