# --------------------------------------# -------------------------- # Math 1225 Laboratory Assignment #4 February 24, 1999 # Treibergs Due March 10, 1999 # ------------------------------------------------------------------ # Taylor Polynomials and approximation of functions. # The Taylor polynomial computed for a function f(x) of order n about the center point a is # given by formula p474. > with(student): > ex1:=taylor(f(x),x=a,6); (2) 2 ex1 := f(a) + D(f)(a) (x - a) + 1/2 D (f)(a) (x - a) (3) 3 (4) 4 + 1/6 D (f)(a) (x - a) + 1/24 D (f)(a) (x - a) (5) 5 6 + 1/120 D (f)(a) (x - a) + O((x - a) ) # The MAPLE instruction to print out the approximating polynomial is "taylor(f(x),x=a,n+1)" # where "f(x)" is a function to be approximated, "x=a" is the point at which the # approximation occurs, and "n" is the order of the polynomial as in the book. "taylor" also # prints the magnitude of the error the polynomial makes near the point "a". For example # "O((x-2)^m)" means that for some constant "k", the function "f(x)" and the polynomial # "p(x)" will differ by an amount "|f(x)-p(x)| < k|x-a|^ m" as "x" tends to "a". # To cut off the "O" part of the expression, we convert it to a polynomial. > convert(ex1,polynom); (2) 2 (3) 3 f(a) + D(f)(a) (x - a) + 1/2 D (f)(a) (x - a) + 1/6 D (f)(a) (x - a) (4) 4 (5) 5 + 1/24 D (f)(a) (x - a) + 1/120 D (f)(a) (x - a) # # For example "taylor(x^2-5*x+3,x=4,2)" is the straight line approximation = the tangent # line to the curve "y=x^2-5*x+3" at "x=4". Lets see: "dy/dx = 2*x-5" so at "x=4" the # slope is "dy/dx=3" and "y=-1". Thus the point slope form of the tangent line is "y+1=3*(x- # 4)" or "y=-1+3*(x-4)". Check: > taylor(x^2-5*x+3,x=4,2); 2 - 1 + 3 (x - 4) + O((x - 4) ) # Lets try it for a specific function at a=.5. To begin with, lets compute the derivatives one by # one. (This section is inspired by Robert Lopez, MAPLE via Calculus, Birkhouser 1994, p.44.) # Then we work out the coefficients one by one. Then assemble the polynomials. > f:=arctan(5*x^2-2*x-1);a:=1/2; 2 f := arctan(5 x - 2 x - 1) a := 1/2 # Then we compute each derivative successively, adding up the Taylor polynomials. > der.0:=f;c.0:=subs(x=a,f);p.0:=c.0; 2 der0 := arctan(5 x - 2 x - 1) c0 := arctan(-3/4) p0 := - arctan(3/4) > der.1:=diff(der.0,x);c.1:=subs(x=a,der.1);p.1:=p.0+c.1*(x-a); 10 x - 2 der1 := --------------------- 2 2 1 + (5 x - 2 x - 1) 48 c1 := ---- 25 48 24 p1 := - arctan(3/4) + ---- x - ---- 25 25 > der.2:=diff(der.1,x);c.2:=subs(x=a,der.2)/2;p.2:=p.1+c.1*(x-a)^2; 2 2 10 (10 x - 2) (5 x - 2 x - 1) der2 := --------------------- - 2 ---------------------------- 2 2 2 2 2 1 + (5 x - 2 x - 1) (1 + (5 x - 2 x - 1) ) 3728 c2 := ---- 625 48 24 48 2 p2 := - arctan(3/4) + ---- x - ---- + ---- (x - 1/2) 25 25 25 > der.3:=diff(der.2,x);c.3:=subs(x=a,der.3)/3!;p.3:=p.2+c.2*(x-a)^3; 2 3 2 2 (5 x - 2 x - 1) (10 x - 2) (10 x - 2) (5 x - 2 x - 1) der3 := - 60 --------------------------- + 8 ----------------------------- 2 2 2 2 2 3 (1 + (5 x - 2 x - 1) ) (1 + (5 x - 2 x - 1) ) 3 (10 x - 2) - 2 ------------------------ 2 2 2 (1 + (5 x - 2 x - 1) ) 169344 c3 := ------ 15625 48 24 48 2 3728 3 p3 := - arctan(3/4) + ---- x - ---- + ---- (x - 1/2) + ---- (x - 1/2) 25 25 25 625 > plot({f,p.(0..3)},x=-1..1,y=-2..2); # Of course this can be built into a loop. Note that f'(1)=0 so that there is new term p.1. > f:=x/(1+x^2);a:=1;der.0:=f;c.0:=subs(x=a,f);p.0:=c.0; x f := ------ 2 1 + x a := 1 x der0 := ------ 2 1 + x c0 := 1/2 p0 := 1/2 > for k from 1 to 3 do;der.k:=diff(der.(k-1),x); c.k:=subs(x=a,der.k)/k!;p.k:=p.(k-1)+c.k*(x-a)^k;od; 2 1 x der1 := ------ - 2 --------- 2 2 2 1 + x (1 + x ) c1 := 0 p1 := 1/2 3 x x der2 := - 6 --------- + 8 --------- 2 2 2 3 (1 + x ) (1 + x ) c2 := -1/4 2 p2 := 1/2 - 1/4 (x - 1) 2 4 x 6 x der3 := 48 --------- - --------- - 48 --------- 2 3 2 2 2 4 (1 + x ) (1 + x ) (1 + x ) c3 := 1/4 2 3 p3 := 1/2 - 1/4 (x - 1) + 1/4 (x - 1) > plot({f,p.(0..3)},x=-3..3,y=-2..2); # Now we let MAPLE do all of the work, > f:=x-> cos(x); f := cos # If the Taylor polynomials are evaluated at "x=0" then they are called Maclaurin polynomials. # Make a set consisiting of the function, and several Maclaurin polynomials. Because cosine is # even, so are the Maclaurin polynomials. This is because the odd derivatives of "f" vanish at # "x=0". Thus "taylor(cos(x),x=0,5)" and "taylor(cos(x),x=0,6)" return the same # polynomial part. Then plot the functions. Figure out which graph is "cos(x)" and which are the # polynomials of degree "m". Note that some coefficients are zero so that p.(2*k)=p.(2*k+1). > S:={f(x),seq(convert(taylor(f(x),x=0,m),polynom),m=1..18)}; 2 4 6 2 4 S := {1, cos(x), 1 - 1/2 x + 1/24 x - 1/720 x , 1 - 1/2 x + 1/24 x , 2 2 4 6 8 10 1 - 1/2 x , 1 - 1/2 x + 1/24 x - 1/720 x + 1/40320 x - 1/3628800 x 12 14 + 1/479001600 x - 1/87178291200 x , 2 4 6 8 10 1 - 1/2 x + 1/24 x - 1/720 x + 1/40320 x - 1/3628800 x 12 + 1/479001600 x , 2 4 6 8 10 1 - 1/2 x + 1/24 x - 1/720 x + 1/40320 x - 1/3628800 x , 2 4 6 8 1 - 1/2 x + 1/24 x - 1/720 x + 1/40320 x , 2 4 6 8 10 1 - 1/2 x + 1/24 x - 1/720 x + 1/40320 x - 1/3628800 x 12 14 16 + 1/479001600 x - 1/87178291200 x + 1/20922789888000 x } > plot(S,x=0..14,y=-7..7,title=`Maclaurin approximations to y = cos x`); # Notice how the higher degree polynomials do a better job wrapping around the cosine function # and therefore do a better job of approximating the function. For the next experiment, consider # insteaad the Taylor polynomials at x=Pi/4. This time, none of the derivatives of cosine vanish # at 3*PI/4 so there are Taylor polynomials of every degree. Notice how well the Taylor # polynomials follow the function this time. MAPLE recognizes "Pi" as the constant. > S:={f(x),seq(convert(taylor(f(x),x=3*Pi/4,m),polynom),m=1..6)}; 1/2 1/2 S := {- 1/2 2 - 1/2 2 (x - 3/4 Pi) + 1/4 %1, cos(x), 1/2 1/2 1/2 - 1/2 2 - 1/2 2 (x - 3/4 Pi), - 1/2 2 , 1/2 1/2 1/2 3 - 1/2 2 - 1/2 2 (x - 3/4 Pi) + 1/4 %1 + 1/12 2 (x - 3/4 Pi) , 1/2 1/2 1/2 3 - 1/2 2 - 1/2 2 (x - 3/4 Pi) + 1/4 %1 + 1/12 2 (x - 3/4 Pi) 1/2 4 - 1/48 2 (x - 3/4 Pi) , 1/2 1/2 1/2 3 - 1/2 2 - 1/2 2 (x - 3/4 Pi) + 1/4 %1 + 1/12 2 (x - 3/4 Pi) 1/2 4 1/2 5 - 1/48 2 (x - 3/4 Pi) - 1/240 2 (x - 3/4 Pi) } 1/2 2 %1 := 2 (x - 3/4 Pi) > plot(S,x=0..14,y=-7..7,title=`Taylor approximations to y = cos x at x=3*Pi/4`); # We try it with a third function y=ln(x). This time, the polynomials do a much worse job of # following the function . > f:=x-> ln(x); f := ln > S:={f(x),seq(convert(taylor(f(x),x=1,m),polynom),m=1..8)}; 2 3 4 S := {0, x - 1 - 1/2 (x - 1) + 1/3 (x - 1) - 1/4 (x - 1) , 2 x - 1 - 1/2 (x - 1) , 2 3 4 5 x - 1 - 1/2 (x - 1) + 1/3 (x - 1) - 1/4 (x - 1) + 1/5 (x - 1) , 2 3 4 5 x - 1 - 1/2 (x - 1) + 1/3 (x - 1) - 1/4 (x - 1) + 1/5 (x - 1) 6 7 - 1/6 (x - 1) + 1/7 (x - 1) , 2 3 4 5 x - 1 - 1/2 (x - 1) + 1/3 (x - 1) - 1/4 (x - 1) + 1/5 (x - 1) 6 - 1/6 (x - 1) , x - 1, 2 3 ln(x), x - 1 - 1/2 (x - 1) + 1/3 (x - 1) } > plot(S,x=0..14,y=-7..7,title=`Taylor approximations to y = ln(x) at x=1`); # The Taylor's remainder Formula (p479) gives the theoretical error made by approximating # the function by its Taylor polynomial p[n] of order n. It estimates |f(b)-p[n](b)| in terms # of a,b, and n. The drawback to the formula is that a bound for "D^(n+1)f(c)" has to be found # for all c between a and b. In case of f(x)=cos(x) this is easy since |D^(n+1)(c)|<= 1 for # all c. Using this estimate, we can compare the theoretical error and the actual error. We # center the polynomials at a and compute f(b)-p[n](b) and the bound for the remainder term # R[n](b). > f:=x->cos(x);a:=2.4;b:=3.0; f := cos a := 2.4 b := 3.0 # Now we make a table showing n, the actual error f(b)-p[n](b) and the estimate for the error # |R[n](c)|<= |b-a|^(n+1)/(n+1)!. Note that the actual error is smaller for all n until the # roundoff errors become significant. > k:='k';DerBound:=1;RBound:=DerBound*abs(a-b)^k/k!; k := k DerBound := 1 k .6 RBound := --- k! > for k from 1 to 6 do; k,f(b)-subs(x=b,convert(taylor(f(x),x=a,k),polynom)),RBound;od; 1, -.2525987811, .6 2, .1526791284, .1800000000 3, .0199482594, .03600000001 4, -.0043684151, .005400000000 5, -.0003864890, .0006480000000 6, .0000512111, .00006480000001 # Lets try to compute the errors for the approximation of ln(x) at b=2. D ln(x) = x^(-1), # and D^k ln(x) = (-1)^(k-1)(k-1)! x^(-k). Thus |D^k ln(c)| <= (k-1)! in the worst # case for 1 k:='k';f:= x-> ln(x); a:= 1.0; b:= 2.0; k := k f := ln a := 1.0 b := 2.0 > DerBound:=(k-1)!;RBound:=DerBound*abs(a-b)^k/k!; DerBound := (k - 1)! k (k - 1)! 1.0 RBound := ------------- k! > `n actual error theoretical error bound`; for k from 2 to 7 do; k,` `,f(b)-subs(x=b,convert(taylor(f(x),x=a,k),polynom)), ` `,RBound; od; n actual error theoretical error bound 2, , -.3068528194, , .5000000000 3, , .1931471806, , .3333333334 4, , -.1401861527, , .2500000000 5, , .1098138473, , .2000000000 6, , -.0901861527, , .1666666667 7, , .0764805140, , .1428571428 > plot({f(x),seq(convert(taylor(f(x),x=a,k),polynom),k=1..6)},x=0..14,y=-7..7,title=`Taylor approximation to y=ln(x) at x=a`); # Trapezoidal rule. # To approximate the area under a curve y=g(x), the computer can subdivide the domain into n # equally spaced intervals and add up the areas of the little trapezoids. For example, if a and a+h # are x-coordinates then the linear function linear=A+ B*x which satisfies l(a)=g(a) and # l(a+h) = g(a+h) can be found by solving for the coefficients A and B. > a:='a';b:='b';f:='f';x:='x'; a := a b := b f := f x := x > linear := A.0 + A.1*x; linear := A0 + A1 x > eq1:= subs(x=a,linear)=f(a); eq1 := A0 + A1 a = f(a) > eq2:=subs(x=a+h,linear)=f(a+h); eq2 := A0 + A1 (a + h) = f(a + h) > solve({eq1,eq2},{A.0,A.1}); f(a) - f(a + h) a f(a) - a f(a + h) + f(a) h {A1 = - ---------------, A0 = ----------------------------} h h > subs(",linear); a f(a) - a f(a + h) + f(a) h (f(a) - f(a + h)) x ---------------------------- - ------------------- h h > int(",x=a..a+h); (a + h) (a f(a) - a f(a + h) + f(a) h + f(a + h) h) 1/2 --------------------------------------------------- h a (a f(a) - a f(a + h) + 2 f(a) h) - 1/2 ---------------------------------- h > T:=simplify("); T := 1/2 h (f(a) + f(a + h)) > sum(subs({a=L+j*(R-L)/n,h=(R-L)/n},T),j=0..n); n / j (R - L) j (R - L) R - L \ ----- (R - L) |f(L + ---------) + f(L + --------- + -----)| \ \ n n n / ) 1/2 ----------------------------------------------------- / n ----- j = 0 # We can compare with the canned trapezoid rule. Our sum agrees with the trapezoid formula # because all of the interior points L+j*(R-L)/n for j=1..n-1 occur in our sum twice! > trapezoid(f(x),x=L..R,n); / /n - 1 \ \ | |----- | | | | \ i (R - L) | | (R - L) |f(L) + 2 | ) f(L + ---------)| + f(R)| | | / n | | | |----- | | \ \i = 1 / / 1/2 -------------------------------------------------- n # There is another formula, the midpoint rule. For each rectangle from x=a to x= a+h one # simply assumes that the height at the midpoint is a good guess for the of the average of the # function f: ie., the average int(f(x),x=a..a+h)/h is approximated by f(a+h/2). Thus # summing, > f:='f';sum(subs({h=(L-R)/n,a=L+j*(L-R)/n},f(a+h/2)*h),j=0..n-1); f := f n - 1 j (L - R) L - R ----- f(L + --------- + 1/2 -----) (L - R) \ n n ) ------------------------------------ / n ----- j = 0 > middlesum(f(x),x=L..R,n); /n - 1 \ |----- | | \ (i + 1/2) (R - L) | (R - L) | ) f(L + -----------------)| | / n | |----- | \i = 0 / ---------------------------------------- n # Finally, lets see how well the various quadrature (numerical integration) formulae work. # Choose a function whose integral we know. int(4/(1+x^2),x=0..1)=Pi. > f:=x->4/(1+x^2); 4 f := x -> ------ 2 1 + x > int(f(x),x=0..1);leftbox(f(x),x=0..1,10); Pi > ` No. Left sum Right Sum Trapezoid Middle Sum Error Midpt`;for mm from 2 to 8 do; mm,evalf(leftsum(f(x),x=0..1,mm)), evalf(rightsum(f(x),x=0..1,mm)),evalf(trapezoid(f(x),x=0..1,mm)), evalf(middlesum(f(x),x=0..1,mm)),evalf(middlesum(f(x),x=0..1,mm)-Pi);od; No. Left sum Right Sum Trapezoid Middle Sum Error Midpt 2, 3.600000000, 2.600000000, 3.100000000, 3.162352941, .020760287 3, 3.456410256, 2.789743589, 3.123076923, 3.150849210, .009256556 4, 3.381176470, 2.881176470, 3.131176471, 3.146800518, .005207864 5, 3.334926114, 2.934926114, 3.134926114, 3.144925864, .003333210 6, 3.303629734, 2.970296401, 3.136963067, 3.143907427, .002314773 7, 3.281048454, 2.995334168, 3.138191311, 3.143293318, .001700664 8, 3.263988495, 3.013988495, 3.138988495, 3.142894730, .001302076 # PROBLEMS FOR LAB ASSIGNMENT 4. # 1. Let f(x)=sin(x). a.) Work out the Taylor polynomial about a=0 (the Maclaurin # polynomial) for f up to third order by hand. b.) Use MAPLE to find the Taylor polynomials for # f up to 17th order. c.) Plot f and p.0 through p.17. d.) Choose a nonzero number b between # -1 and 1. Work out the error made using the remainder formula. Make a table showing k, # the actual error p.k(b)-f(b), and the theoretically computed error for k=0..17. # 2. Let f(x)=exp(x). a.) Work out the Taylor polynomial about a=1 for f up to third order by # hand. b.) Use MAPLE to find the Taylor polynomials for f up to 9th order. c.) Plot f and p.0 # through p.9. d.) Choose number b strictly between 1 and 2. Work out the error made using # the remainder formula. Make a table showing k, the actual error |p.k(b)-f(b)|, and the # theoretically computed error for k=0..9. # 3. Given three points a, a+h and a+2h on the x-axis, and a function f(x), use MAPLE to # find the coefficients of a quadratic function q(x) =Q.0 + Q.1*x + Q.2*x^2 such that # q(a)=f(a), q(a+h)=f(a+h) and q(a+2*h)=f(a+2*h). This f(x) is pretty well approximated # by q(x) if the third derivativatives are not large. Choose a function f(x) and numbers a, # and small h. Plot your f and the corresponding q(x) over the interval a-h..a+3*h. Use # MAPLE to find the integral "int(q(x),x=a..a+2*h);" This approximation of the integral # "int(f(x),x=a..a+2*h)" called Simpson's rule. Divide the interval [L,R] into n equal # subintervals. For the j-th subinterval, set a= L+j*(R-L)/n and h=(L-R)/(2*n) and add # up the resulting areas for j=0..n-1. Compare your answer to the result of # "simplson(f(x),x=L..R,n);". # 4. Let f(x)=sqrt(4-x^2); Find int(f(x),x=0..2); Make a table comparing the trapezoid # rule, midpoint rule and Simpson's rule for n=2,4,6,8,...,1024 rectangles approximating the # integral of this function. Print out the integrals and errors made by these three rules. Repeat # with another function f(x) and interval [L,R] of your choice. Compare the three quadrature # rules in how well they of approximate the integral.