# 1.33 Lesson 33

```# 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.

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

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 **

** 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.

(.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.

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 )