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