>   

>    restart;

>   

Adding Numbers of Different Sizes

>    a := 10^8;

a := 100000000

>    b := 0.4;

b := .4

>    Digits := 8;     # good values are 10, 9, 8, 7, 6

Digits := 8

>    s := a:

>    for i from 1 while s=a do

>      bb := add( b, k=1..i );

>      s := a+bb;

>      printf("%6d  %6.2f   %12.2f\n", i, bb, s );

>    end do: 

     1     .40   100000000.00

     2     .80   100000000.00

     3    1.20   100000000.00

     4    1.60   100000000.00

     5    2.00   100000000.00

     6    2.40   100000000.00

     7    2.80   100000000.00

     8    3.20   100000000.00

     9    3.60   100000000.00

    10    4.00   100000000.00

    11    4.40   100000000.00

    12    4.80   100000000.00

    13    5.20   100000010.00

>    Digits := 10;    # restore to default value

Digits := 10

>   

Subtracting Numbers that are Very Close

>    f := x -> x^4;

f := proc (x) options operator, arrow; x^4 end proc

>    a := 2;

a := 2

>    Digits := 20;     # good values are integers up to 20

Digits := 20

>    h := 1:

>    for n from 0 to 20 do

>      h := 10^(-n);

>      q := fnormal(f(2+h)-f(2))/h;

>      printf("%3d %25.20f %10.5f %25.20f %25.20f\n", n, h, f(2), f(2+h), q );

>    end do: 

  0    1.00000000000000000000   16.00000   81.00000000000000000000   65.00000000000000000000

  1     .10000000000000000000   16.00000   19.44810000000000000000   34.48100000000000000000

  2     .01000000000000000000   16.00000   16.32240801000000000000   32.24080100000000000000

  3     .00100000000000000000   16.00000   16.03202400800100000000   32.02400800100000000000

  4     .00010000000000000000   16.00000   16.00320024000800010000   32.00240008000100000000

  5     .00001000000000000000   16.00000   16.00032000240000800000   32.00024000080000100000

  6     .00000100000000000000   16.00000   16.00003200002400000800   32.00002400000800000100

  7     .00000010000000000000   16.00000   16.00000320000024000000   32.00000240000008000000

  8     .00000001000000000000   16.00000   16.00000032000000240000   32.00000024000000080000

  9     .00000000100000000000   16.00000   16.00000003200000002400   32.00000002400000000800

 10     .00000000010000000000   16.00000   16.00000000320000000000   32.00000000240000000000

 11     .00000000001000000000   16.00000   16.00000000032000000000   32.00000000024000000000

 12     .00000000000100000000   16.00000   16.00000000003200000000   32.00000000002400000000

 13     .00000000000010000000   16.00000   16.00000000000320000000   32.00000000000240000000

 14     .00000000000001000000   16.00000   16.00000000000032000000   32.00000000000024000000

 15     .00000000000000100000   16.00000   16.00000000000003200000   32.00000000000002400000

 16     .00000000000000010000   16.00000   16.00000000000000320000   32.00000000000000240000

 17     .00000000000000001000   16.00000   16.00000000000000032000   32.00000000000000024000

 18     .00000000000000000100   16.00000   16.00000000000000003200   32.00000000000000002400

 19     .00000000000000000010   16.00000   16.00000000000000000300   32.00000000000000000200

 20     .00000000000000000001   16.00000   16.00000000000000000000    0.00000000000000000000

>    Digits := 10;     # restore to default value

Digits := 10

>   

Approximation by Taylor/Maclaurin Polynomials

>    Digits := 10;

Digits := 10

# 4

Enter information about problem

>    f := x -> tan(x);   # function to be approximated

f := tan

>    X := 0.12;          # point where function is to be evaluated

X := .12

>    a := 0;             # expansion point (Maclaurin -> a=0)

a := 0

>    n := 4;             # order of expansion

n := 4

>   

Find Taylor/Maclaurin polynomial and error term

>    p[n] := convert( taylor( f(x), x=a, n+1 ), polynom ):

>    P[n] := unapply( p[n], x );  # Taylor polynomial

P[4] := proc (x) options operator, arrow; x+1/3*x^3 end proc

>    m[n] := simplify((D@@(n+1))(f)(c)):

>    r[n] := unapply( m[n]*(x-a)^(n+1)/(n+1)!, (c,x) );   # exact remainder

r[4] := proc (c, x) options operator, arrow; 1/15*(1+tan(c)^2)*(15*tan(c)^4+15*tan(c)^2+2)*x^5 end proc

>   

Find bound for error term

>    plot( abs(m[n]), c=min(a,X)..max(a,X) );   # looking for maximum value between X and a

[Maple Plot]

>    M[n] := 20;  # bound for coefficient

M[4] := 20

>    R[n] := unapply( M[n]*(x-a)^(n+1)/(n+1)!, x );

R[4] := proc (x) options operator, arrow; 1/6*x^5 end proc

>    plot( [-R[n](x), f(x)-P[n](x), R[n](x)], x=0..0.12,

>          color=[pink,blue,pink],

>          legend=["error bound","actual difference","error bound"] );

[Maple Plot]

>   

Compare approximate value, exact value, and approximation error

>    evalf( [P[n](X), f(X), P[n](X)-f(X)] );

[.1205760000, .1205793372, -.33372e-5]

>   

# 14

Enter information about problem

>    f := x -> sqrt(x);   # function to be approximated

f := sqrt

>    X := 5;          # point where function is to be evaluated

X := 5

>    a := 4;             # expansion point (Maclaurin -> a=0)

a := 4

>    n := 4;             # order of expansion

n := 4

>   

Find Taylor/Maclaurin polynomial and error term

>    p[n] := convert( taylor( f(x), x=a, n+1 ), polynom ):

>    P[n] := unapply( p[n], x );  # Taylor polynomial

P[4] := proc (x) options operator, arrow; 1+1/4*x-1/64*(x-4)^2+1/512*(x-4)^3-5/16384*(x-4)^4 end proc

>    m[n] := simplify((D@@(n+1))(f)(c)):

>    r[n] := unapply( m[n]*(x-a)^(n+1)/(n+1)!, (c,x) );   # exact remainder

r[4] := proc (c, x) options operator, arrow; 7/256*1/c^(9/2)*(x-4)^5 end proc

>   

Find bound for error term

>    plot( abs(m[n]), c=min(a,X)..max(a,X) );   # looking for maximum value between X and a

[Maple Plot]

>    M[n] := 0.15;  # bound for coefficient

M[4] := .15

>    R[n] := unapply( M[n]*(x-a)^(n+1)/(n+1)!, x );

R[4] := proc (x) options operator, arrow; .1250000000e-2*(x-4)^5 end proc

>    plot( [-R[n](x), f(x)-P[n](x), R[n](x)], x=min(a,X)..max(a,X),

>          color=[pink,blue,pink],

>          legend=["error bound","actual difference","error bound"] );

[Maple Plot]

>   

Compare approximate value, exact value, and approximation error

>    evalf( [P[n](X), f(X), P[n](X)-f(X)] );

[2.236022949, 2.236067977, -.45028e-4]

>   

# 38

Enter information about problem

>    f := x -> sin(x);   # function to be approximated

f := sin

>    X := Pi/8;          # point where function is to be evaluated

X := 1/8*Pi

>    a := Pi/4;             # expansion point (Maclaurin -> a=0)

a := 1/4*Pi

>    n := 3;             # order of expansion

n := 3

>   

Find Taylor/Maclaurin polynomial and error term

>    p[n] := convert( taylor( f(x), x=a, n+1 ), polynom ):

>    P[n] := unapply( p[n], x );  # Taylor polynomial

P[3] := proc (x) options operator, arrow; 1/2*2^(1/2)+1/2*2^(1/2)*(x-1/4*Pi)-1/4*2^(1/2)*(x-1/4*Pi)^2-1/12*2^(1/2)*(x-1/4*Pi)^3 end proc

>    m[n] := simplify((D@@(n+1))(f)(c)):

>    r[n] := unapply( m[n]*(x-a)^(n+1)/(n+1)!, (c,x) );   # exact remainder

r[3] := proc (c, x) options operator, arrow; 1/24*sin(c)*(x-1/4*Pi)^4 end proc

>   

Find bound for error term

>    plot( abs(m[n]), c=min(a,X)..max(a,X) );   # looking for maximum value between X and a

[Maple Plot]

>    M[n] := sqrt(2)/2;  # bound for coefficient

M[3] := 1/2*2^(1/2)

>    R[n] := unapply( M[n]*(x-a)^(n+1)/(n+1)!, x );

R[3] := proc (x) options operator, arrow; 1/48*2^(1/2)*(x-1/4*Pi)^4 end proc

>    plot( [-R[n](x), f(x)-P[n](x), R[n](x)], x=min(a,X)..max(a,X),

>          color=[pink,blue,pink],

>          legend=["error bound","actual difference","error bound"] );

[Maple Plot]

>   

Compare approximate value, exact value, and approximation error

>    evalf( [P[n](X), f(X), P[n](X)-f(X)] );

[.3820411832, .3826834325, -.6422493e-3]

>   

# 42

Enter information about problem

>    f := x -> ln((1+x)/(1-x));   # function to be approximated

f := proc (x) options operator, arrow; ln((1+x)/(1-x)) end proc

>    #X := 5;          # point where function is to be evaluated

>    a := 0;             # expansion point (Maclaurin -> a=0)

a := 0

>    n := 4;             # order of expansion

n := 4

>   

Find Taylor/Maclaurin polynomial and error term

>    p[n] := convert( taylor( f(x), x=a, n+1 ), polynom ):

>    P[n] := unapply( p[n], x );  # Taylor polynomial

P[4] := proc (x) options operator, arrow; 2*x+2/3*x^3 end proc

>    m[n] := simplify((D@@(n+1))(f)(c)):

>    r[n] := unapply( m[n]*(x-a)^(n+1)/(n+1)!, (c,x) );   # exact remainder

r[4] := proc (c, x) options operator, arrow; -2/5*(1+10*c^2+5*c^4)/(1+c)^5/(-1+c)^5*x^5 end proc

>   

Find bound for error term

>    plot( abs(m[n]), c=-0.5..0.5 );   # looking for maximum value between X and a

[Maple Plot]

>    M[n] := 800;  # bound for coefficient

M[4] := 800

>    R[n] := unapply( M[n]*(x-a)^(n+1)/(n+1)!, x );

R[4] := proc (x) options operator, arrow; 20/3*x^5 end proc

>    plot( [-R[n](x), f(x)-P[n](x), R[n](x)], x=-0.5..0.5,

>          color=[pink,blue,pink],

>          legend=["error bound","actual difference","error bound"] );

[Maple Plot]

>   

>   

>