# 3.3 Assignment 3 Solutions(Israel's)

```# Math 210 Sec. 201
# Solutions to Assignment 3
#
# 1.
> r:= cos(theta)^2:
> L:= Int(sqrt(r^2 + diff(r,theta)^2), theta = -Pi/2 .. Pi/2);

1/2 Pi
/
|                4               2           2 1/2
L :=     |     (cos(theta)  + 4 cos(theta)  sin(theta) )    dtheta
|
/
- 1/2 Pi
> L:= int(cos(theta)*sqrt(cos(theta)^2+4*sin(theta)^2), theta = -Pi/2 .. Pi/2);

1/2 Pi
/
|                           2               2 1/2
L :=     |     cos(theta) (cos(theta)  + 4 sin(theta) )    dtheta
|
/
- 1/2 Pi
> with(student):

> changevar(x=sin(theta),L,x);

1/2          1/2
2 + 1/3 3    arcsinh(3   )
> convert(",ln);

1/2         1/2
2 + 1/3 3    ln(2 + 3   )
> evalf(");

2.760345997
# 2. (a)
> f:= x -> x^11:

# The following definitions are copied from Lesson 18.
> h:= n ->(b-a)/n:
> X:= (i,n) -> a + i*(b-a)/n:
> midrule:= n ->
>    h(n)*sum(f((X(i,n)+X(i+1,n))/2), i = 0 .. n-1);
> traprule:= n ->
>    h(n) * sum((f(X(i,n))+f(X(i+1,n)))/2, i=0..n-1);
> simpson:= n ->
>    h(n) * sum(1/3*f(X(2*i,n))+4/3*f(X(2*i+1,n))+1/3*f(X(2*i+2,n)), i=0..n/2-1);
> J:= int(f(x), x=a..b);

/n - 1                                 \
|-----                                 |
| \                                    |
midrule := n -> h(n) |  )   f(1/2 X(i, n) + 1/2 X(i + 1, n))|
| /                                    |
|-----                                 |
\i = 0                                 /

/n - 1                                      \
|-----                                      |
| \                                         |
traprule := n -> h(n) |  )   (1/2 f(X(i, n)) + 1/2 f(X(i + 1, n)))|
| /                                         |
|-----                                      |
\i = 0                                      /

simpson := n -> h(n)

/1/2 n - 1                                                                 \
|  -----                                                                   |
|   \                                                                      |
|    )     (1/3 f(X(2 i, n)) + 4/3 f(X(2 i + 1, n)) + 1/3 f(X(2 i + 2, n)))|
|   /                                                                      |
|  -----                                                                   |
\  i = 0                                                                   /

b
/
|
J :=  |  f(x) dx
|
/
a
# Now define ME, TE and SE as functions of n.
> ME:= unapply(J - midrule(n), n);

10            8  3
12         12                      a b             a  b
ME := n -> 1/12 b   - 1/12 a   - (b - a) (7665/2048 ----- + 1705/96 -----
9               5
n               n

3  8                              8  3              10
a  b          9  2                a  b              a   b
- 48895/1024 ----- + 1/12 b  a  n + 63875/2048 ----- + 7665/2048 -----
7                                 9                 9
n                                 n                 n

11               4  7                           9  2
a                a  b          7  4             a  b
+ 77/64 --- - 38325/1024 ----- + 1/12 a  b  n - 1705/96 -----
3                 9                              5
n                 n                              n

5  6             7  4                              5  6
a  b             a  b          6  5                a  b
+ 341/192 ----- - 1705/192 ----- + 1/12 a  b  n + 17885/1024 -----
5                5                                 9
n                n                                 n

3  8               10              7  4               2  9
a  b             a b               a  b               a  b
+ 63875/2048 ----- + 1705/192 ----- + 23749/512 ----- + 29337/1024 -----
9                5                 7                  7
n                n                 n                  n

8  3                                          6  5
a  b            10           10               a  b
- 77/64 ----- + 1/12 a b   n + 1/12 a   b n + 341/192 -----
3                                             5
n                                             n

7  4               2  9           10                4  7
a  b               a  b           a   b             a  b
- 38325/1024 ----- - 89425/6144 ----- - 231/64 ----- + 23749/512 -----
9                  9              3                 7
n                  n              n                 n

8  3            9  2          10
9  2                a  b            b  a          a   b
+ 1/12 a  b  n - 48895/1024 ----- - 1705/96 ----- + 11/24 -----
7               5            n
n               n

10             10               6  5           2  9
a b              a   b            a  b           a  b
+ 11/24 ----- + 1705/192 ----- - 9779/512 ----- + 231/64 -----
n                 5                7              3
n                n              n

6  5            8  3          3  8
a  b            b  a          a  b          3  8
+ 17885/1024 ----- + 1705/96 ----- - 77/64 ----- + 1/12 a  b  n
9               5             3
n               n             n

4  7              10                                9  2
a  b              b   a         4  7                a  b
- 1705/192 ----- - 9779/1024 ----- + 1/12 a  b  n - 89425/6144 -----
5                 7                                 9
n                 n                                 n

10                 9  2             5  6
b   a              a  b             a  b          8  3
- 231/64 ----- + 29337/1024 ----- - 9779/512 ----- + 1/12 a  b  n
3                  7                7
n                  n                n

9  2                             10            11
a  b          5  6               a   b         a             11
+ 231/64 ----- + 1/12 a  b  n - 9779/1024 ----- - 11/24 --- + 1/12 n b
3                                7           n
n                                n

11              11              11              11
a               a               a               b
- 341/192 --- + 1397/1024 --- - 2555/6144 --- - 2555/6144 ---
5               7               9               9
n               n               n               n

11            11          11          11
b             b           b           b           11
+ 1397/1024 --- - 341/192 --- + 77/64 --- - 11/24 --- + 1/12 a   n)/n
7             5           3          n
n             n           n
> TE:= unapply(J - traprule(n), n);

10         8  3
12         12                   a b          a  b
TE := n -> 1/12 b   - 1/12 a   - (b - a) (- 15/4 ----- - 55/3 -----
9            5
n            n

3  8                         8  3         10           11
a  b          9  2           a  b         a   b        a
+ 385/8 ----- + 1/12 b  a  n - 125/4 ----- - 15/4 ----- - 11/8 ---
7                            9            9           3
n                            n            n           n

4  7                        9  2         5  6         7  4
a  b          7  4          a  b         a  b         a  b
+ 75/2 ----- + 1/12 a  b  n + 55/3 ----- - 11/6 ----- + 55/6 -----
9                           5            5            5
n                           n            n            n

5  6          3  8           10          7  4
6  5          a  b          a  b         a b           a  b
+ 1/12 a  b  n - 35/2 ----- - 125/4 ----- - 55/6 ----- - 187/4 -----
9             9            5             7
n             n            n             n

2  9         8  3                                       6  5
a  b         a  b            10           10            a  b
- 231/8 ----- + 11/8 ----- + 1/12 a b   n + 1/12 a   b n - 11/6 -----
7            3                                          5
n            n                                          n

7  4           2  9         10            4  7
a  b           a  b         a   b         a  b          9  2
+ 75/2 ----- + 175/12 ----- + 33/8 ----- - 187/4 ----- + 1/12 a  b  n
9              9            3             7
n              n            n             n

8  3         9  2          10              10         10
a  b         b  a          a   b         a b          a   b
+ 385/8 ----- + 55/3 ----- - 11/12 ----- - 11/12 ----- - 55/6 -----
7            5            n             n             5
n            n                                        n

6  5         2  9         6  5         8  3         3  8
a  b         a  b         a  b         b  a         a  b
+ 77/4 ----- - 33/8 ----- - 35/2 ----- - 55/3 ----- + 11/8 -----
7            3            9            5            3
n            n            n            n            n

4  7         10                            9  2
3  8          a  b         b   a         4  7            a  b
+ 1/12 a  b  n + 55/6 ----- + 77/8 ----- + 1/12 a  b  n + 175/12 -----
5            7                             9
n            n                             n

10            9  2         5  6                        9  2
b   a         a  b         a  b          8  3          a  b
+ 33/8 ----- - 231/8 ----- + 77/4 ----- + 1/12 a  b  n - 33/8 -----
3             7            7                           3
n             n            n                           n

10            11                      11         11
5  6          a   b         a             11        a          a
+ 1/12 a  b  n + 77/8 ----- + 11/12 --- + 1/12 n b   + 11/6 --- - 11/8 ---
7           n                        5          7
n                                    n          n

11         11         11         11         11          11
a          b          b          b          b           b
+ 5/12 --- + 5/12 --- - 11/8 --- + 11/6 --- - 11/8 --- + 11/12 ---
9          9          7          5          3          n
n          n          n          n          n

11
+ 1/12 a   n)/n
> SE:= unapply(J - simpson(n), n);

10           8  3
12         12                 a b            a  b
SE := n -> 1/12 b   - 1/12 a   - (b - a) (1275 ----- + 1100/3 -----
9              5
n              n

3  8                         8  3         10           11
a  b          9  2           a  b         a   b        a
- 8085/2 ----- + 1/12 b  a  n + 10625 ----- + 1275 ----- + 11/2 ---
7                            9            9           3
n                            n            n           n

4  7                          9  2          5  6          7  4
a  b          7  4            a  b          a  b          a  b
- 12750 ----- + 1/12 a  b  n - 1100/3 ----- + 110/3 ----- - 550/3 -----
9                             5             5             5
n                             n             n             n

5  6          3  8            10         7  4
6  5          a  b          a  b          a b          a  b
+ 1/12 a  b  n + 5950 ----- + 10625 ----- + 550/3 ----- + 3927 -----
9             9             5            7
n             n             n            n

2  9         8  3                                        6  5
a  b         a  b            10           10             a  b
+ 4851/2 ----- - 11/2 ----- + 1/12 a b   n + 1/12 a   b n + 110/3 -----
7            3                                           5
n            n                                           n

7  4            2  9         10           4  7
a  b            a  b         a   b        a  b          9  2
- 12750 ----- - 14875/3 ----- - 33/2 ----- + 3927 ----- + 1/12 a  b  n
9               9            3            7
n               n            n            n

8  3           9  2          10           6  5         2  9
a  b           b  a          a   b        a  b         a  b
- 8085/2 ----- - 1100/3 ----- + 550/3 ----- - 1617 ----- + 33/2 -----
7              5             5            7            3
n              n             n            n            n

6  5           8  3         3  8                         4  7
a  b           b  a         a  b          3  8           a  b
+ 5950 ----- + 1100/3 ----- - 11/2 ----- + 1/12 a  b  n - 550/3 -----
9              5            3                            5
n              n            n                            n

10                             9  2         10             9  2
b   a         4  7             a  b         b   a          a  b
- 1617/2 ----- + 1/12 a  b  n - 14875/3 ----- - 33/2 ----- + 4851/2 -----
7                              9            3              7
n                              n            n              n

5  6                        9  2                          10
a  b          8  3          a  b          5  6            a   b
- 1617 ----- + 1/12 a  b  n + 33/2 ----- + 1/12 a  b  n - 1617/2 -----
7                           3                             7
n                           n                             n

11          11          11          11          11
11         a           a           a           b           b
+ 1/12 n b   - 110/3 --- + 231/2 --- - 425/3 --- - 425/3 --- + 231/2 ---
5           7           9           9           7
n           n           n           n           n

11         11
b          b           11
- 110/3 --- + 11/2 --- + 1/12 a   n)/n
5          3
n          n
# b) To plot these, we need to take particular values of a and b,
# e.g. a=0, b=1.
> a:= 0; b:= 1;
> plot({n^2 * ME(n), n^2* TE(n)}, n = 1 .. 20);

a := 0

b := 1
# The Midpoint error is on the top, and the Trapezoid error is on the bottom.  We see that the Midpoint
# error has about half the absolute value of the Trapezoid error, and the opposite sign.
# The Midpoint error is approximately 0.45/n^2, and the Trapezoid error approximately -0.9/n^2.

> plot(n^4*SE(n), n = 2 .. 40);

# The Simpson's Rule error is approximately -5.4/n^4 for large n.
#
# (c) For this we need a and b to be symbolic again.
> a:= 'a': b:= 'b':

> Mcoeffs:= [seq(coeff(expand(ME(n)),n,-2*jj), jj=1..5)];

Mcoeffs := [

11     11    11   10  2    11   10  2    11   12    11   11      11   12
---- b a   + ---- b   a  - ---- a   b  + ---- b   - ---- b   a - ---- a  ,
12           24            24            24         12           24

77     11    77   8  4    77   9  3    77   11     231  10  2    77   9  3
- ---- b a   + ---- a  b  - ---- a  b  + ---- b   a - --- b   a  + ---- b  a
16           64           16           16           32           16

231  10  2    77   12    77   12    77   8  4
+ --- a   b  + ---- a   - ---- b   - ---- b  a ,
32           64         64         64

1705  8  4   1705  10  2   1705  9  3   341  7  5   1705  8  4   1705  9  3
---- b  a  + ---- b   a  - ---- b  a  + --- a  b  - ---- a  b  + ---- a  b
64           64            48           32          64           48

1705  10  2   341  11     341  5  7   341  12   341  12   341    11
- ---- a   b  - --- b   a - --- a  b  + --- b   - --- a   + --- b a  ,
64            32          32         192       192        32

1397    11   96393  8  4   9779  9  3   1397  11     9779  9  3   96393  8  4
- ---- b a   + ----- a  b  - ---- a  b  + ---- b   a + ---- b  a  - ----- b  a
128          1024          128          128          128          1024

9779  10  2   9779  10  2   4191  5  7   4191  7  5   1397  12   1397  12
+ ---- a   b  - ---- b   a  + ---- a  b  - ---- a  b  - ---- b   + ---- a
256           256           64           64          1024       1024

,

12775    11   12775  11     140525  8  4   28105  5  7   140525  9  3
----- b a   - ----- b   a + ------ b  a  - ----- a  b  - ------ b  a
3072          3072          2048           512           3072

28105  10  2   28105  7  5   140525  8  4   140525  9  3   28105  10  2
+ ----- b   a  + ----- a  b  - ------ a  b  + ------ a  b  - ----- a   b
1536           512           2048           3072           1536

2555  12   2555  12
+ ---- b   - ---- a                                                      ]
6144       6144
> shouldbe := [ seq((b-a)^(2*jj)*((D@@(2*jj-1))(f)(b) - (D@@(2*jj-1))(f)(a)), jj = 1 .. 5) ];

2      10       10          4       8        8
shouldbe := [(b - a)  (11 b   - 11 a  ), (b - a)  (990 b  - 990 a ),

6         6          6          8           4            4
(b - a)  (55440 b  - 55440 a ), (b - a)  (1663200 b  - 1663200 a ),

10            2             2
(b - a)   (19958400 b  - 19958400 a )]
> Mcs:= seq(normal(Mcoeffs[kk]/shouldbe[kk]), kk = 1 .. 5);

31         127         73
Mcs := 1/24, -7/5760, ------, - ---------, ----------
967680    154828800  3503554560
> Tcoeffs:= [seq(coeff(expand(TE(n)),n,-2*jj), jj=1..5)];

Tcoeffs := [

11      11   12    11   10  2    11   12    11   10  2           11
11/6 b   a - ---- b   + ---- a   b  + ---- a   - ---- b   a  - 11/6 b a  ,
12         12            12         12

12         9  3         10  2         8  4         9  3         11
11/8 b   + 11/2 a  b  - 33/4 a   b  + 11/8 b  a  - 11/2 b  a  - 11/2 b   a

10  2         8  4           11         12
+ 33/4 b   a  - 11/8 a  b  + 11/2 b a   - 11/8 a  ,

9  3         10  2         8  4         12         8  4          9  3
- 110/3 a  b  + 55/2 a   b  + 55/2 a  b  + 11/6 a   - 55/2 b  a  + 110/3 b  a

7  5         10  2         12         11       5  7       11
- 11 a  b  - 55/2 b   a  - 11/6 b   - 11 b a   + 11 a  b  + 11 b   a,

12       9  3       5  7         12       9  3         10  2       11
11/8 b   + 77 a  b  - 66 a  b  - 11/8 a   - 77 b  a  - 77/2 a   b  - 11 b   a

8  4          8  4         11       7  5         10  2
+ 759/8 b  a  - 759/8 a  b  + 11 b a   + 66 a  b  + 77/2 b   a ,

8  4          9  3         11           12         12          9  3
- 275/4 b  a  + 275/6 b  a  + 25/6 b   a - 5/12 b   + 5/12 a   - 275/6 a  b

7  5         10  2          8  4         10  2           11
- 55 a  b  - 55/3 b   a  + 275/4 a  b  + 55/3 a   b  - 25/6 b a

5  7
+ 55 a  b                                                              ]
> Tcs:= seq(normal(Tcoeffs[kk]/shouldbe[kk]), kk = 1 .. 5);

Tcs := -1/12, 1/720, -1/30240, 1/1209600, -1/47900160
> Scoeffs:= [seq(coeff(expand(SE(n)),n,-2*jj), jj=1..5)];

Scoeffs := [0,

8  4         11       10  2         8  4       9  3       9  3
- 11/2 b  a  - 22 b a   + 33 a   b  + 11/2 a  b  - 22 a  b  + 22 b  a

11         10  2         12         12
+ 22 b   a - 33 b   a  - 11/2 b   + 11/2 a  ,

10  2        11             9  3        5  7          11        8  4
550 b   a  - 220 b   a + 2200/3 a  b  - 220 a  b  + 220 b a   + 550 b  a

10  2        8  4        7  5           9  3          12          12
- 550 a   b  - 550 a  b  + 220 a  b  - 2200/3 b  a  + 110/3 b   - 110/3 a

,

11         10  2         9  3            8  4         10  2        11
- 924 b a   + 3234 a   b  - 6468 a  b  - 15939/2 b  a  - 3234 b   a  + 924 b   a

8  4         7  5         9  3         5  7          12
+ 15939/2 a  b  - 5544 a  b  + 6468 b  a  + 5544 a  b  - 231/2 b

12
+ 231/2 a  ,

8  4            10  2            9  3          5  7           11
23375 b  a  + 18700/3 b   a  + 46750/3 a  b  - 18700 a  b  - 4250/3 b   a

10  2             11            9  3          7  5          8  4
- 18700/3 a   b  + 4250/3 b a   - 46750/3 b  a  + 18700 a  b  - 23375 a  b

12          12
+ 425/3 b   - 425/3 a

]
> Scs:= seq(normal(Scoeffs[kk]/shouldbe[kk]), kk = 1 .. 5);

17
Scs := 0, -1/180, 1/1512, -1/14400, -------
2395008
# (d)  Let's try it for x^13.
> f:= x -> x^13:

# We must repeat the definitions of J, ME, TE, SE etc.  The coefficients should go up to jj=6 this time.
> J:= int(f(x), x=a..b):
> ME:= unapply(J - midrule(n), n):
> TE:= unapply(J - traprule(n), n):
> SE:= unapply(J - simpson(n), n):
> Mcoeffs:= [seq(coeff(expand(ME(n)),n,-2*jj), jj=1..6)]:
> Tcoeffs:= [seq(coeff(expand(TE(n)),n,-2*jj), jj=1..6)]:
> Scoeffs:= [seq(coeff(expand(SE(n)),n,-2*jj), jj=1..6)]:
> shouldbe := [ seq((b-a)^(2*jj)*((D@@(2*jj-1))(f)(b) - (D@@(2*jj-1))(f)(a)), jj = 1 .. 6) ]:
> Mcs:= seq(normal(Mcoeffs[kk]/shouldbe[kk]), kk = 1 .. 6);
> Tcs:= seq(normal(Tcoeffs[kk]/shouldbe[kk]), kk = 1 .. 6);
> Scs:= seq(normal(Scoeffs[kk]/shouldbe[kk]), kk = 1 .. 6);

31         127         73             1414477
Mcs := 1/24, -7/5760, ------, - ---------, ----------, - ----------------
967680    154828800  3503554560    2678117105664000

691
Tcs := -1/12, 1/720, -1/30240, 1/1209600, -1/47900160, -------------
1307674368000

17         21421
Scs := 0, -1/180, 1/1512, -1/14400, -------, - -----------
2395008    29719872000
# (e)
> a:= 0: b:= 1: f:= exp:
> J:= int(f(x), x=a..b):
> TE:= unapply(J - traprule(n), n);

exp(1) (1 + exp(1/n))       1 + exp(1/n)
1/2 --------------------- - 1/2 ------------
exp(1/n) - 1           exp(1/n) - 1
TE := n -> exp(1) - 1 - --------------------------------------------
n
> formula:= n -> sum( Tcs[j] *(b-a)^(2*j)*((D@@(2*j-1))(f)(b) - (D@@(2*j-1))(f)(a))/ n^(2*j), j = 1 .. 5);

5
-----               (2 j)   (2 j - 1)          (2 j - 1)
\    Tcs[j] (b - a)      (D         (f)(b) - D         (f)(a))
formula := n ->   )   ---------------------------------------------------------
/                               (2 j)
-----                           n
j = 1
> Digits:= 30:
> seq(evalf(nn^12 * (TE(nn)-formula(nn))), nn = 1 .. 10);

Digits := 30

-9                       -9                     -9
.88554587380646207944*10  , .90226022060161953*10  , .905424963218222*10  ,

-9                  -9                 -9
.9065378753292*10  , .907053921534*10  , .90733448936*10  ,

-9                -9               -9              -9
.9075037465*10  , .9076136333*10  , .907688984*10  , .90774289*10
# It certainly appears to be approaching a limit, perhaps somewhere around .908 * 10^(-9).
#
# (f)  Since all derivatives of f(x) = 1/(2 + cos(x)) are the same at 0 and at 2 Pi, the formula would say that
# the error is O(n^(-m)) for every m, i.e. it decreases faster than any power of m.
> n:= 20: a:= 0: b:= 2*Pi: f:= x -> 1/(2+cos(x)):
> J:= int(f(x), x=a .. b);

1/2
J := 2/3 3    Pi
> evalf(J - traprule(n));

-10
-.2640573654902402695*10
> evalf(J - midrule(n));

-10
.2640573654883181638*10
> evalf(J - simpson(n));

-5
.461370711086269911080003*10
# The Trapezoid and Midpoint Rule errors are almost exactly equal in magnitude, and both very small.
# The Simpson's Rule error is much larger.  This is rather surprising at first sight.  The point is that the
# terms that are the main contributions to the errors for most functions happen to be 0 for this f.  The error
# that there is, goes rapidly to 0 as
# n grows.  There is no reason to expect the Midpoint Rule error to be approximately -1/2 times the
# Trapezoid Rule error, as was the case for those terms that happened to be 0 in this case.  So the Simpson's
# Rule error for n = 20, which is a 2/3 : 1/3 mixture of the Midpoint and Trapezoid Rule errors for n=10, is
# of the same order of magnitude as those errors.  But these are much larger than the Midpoint and
# Trapezoid Rule errors for n=20, because the errors converge to 0 so rapidly as n -> infinity.
#
# 3.(a)  I'll keep most of the definitions from question 2, changing what is needed.  Other definitions come
# from Lesson 19.  I'll use Digits = 20, because some of the errors are very small.
> a:= 1: b:= 3: n:= 'n': Digits:= 20:
> f:= x -> 1/(x^2 - 2*x + 5):
> J:= int(f(x), x=a .. b); Je:=evalf(J);

J := 1/8 Pi

Je := .39269908169872415481
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:

> for kk from 2 to 5 do
>   for jj from kk to 5 do
>     R[jj,kk]:= (4^(kk-1)*R[jj,kk-1]-R[jj-1,kk-1])/(4^(kk-1)-1) od od;
> [seq(R[jj,1],jj=1..5)]; [seq(R[jj,2],jj=2..5)];\

> [seq(R[jj,3],jj=3..5)]; [seq(R[jj,4],jj=4..5)];
> R[5,5];

[.38750000000000000000, .39139705882352941176, .39237356181138612617,

.39261770150517361183, .39267873664687180415]

[.39269607843137254900, .39269906280733836430, .39269908140310277370,

.39269908169410453493]

[.39269926176573608532, .39269908264282040100, .39269908171350465235]

[.39269907979959951713, .39269908169875360872]

.39269908170620127181
# (b) For int(x^11, x=0 .. 1):
> a:= 0: b:= 1: f:= x -> x^11:
> J:= int(f(x), x=a..b); Je:= evalf(J);

J := 1/12

Je := .083333333333333333333
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:
> for kk from 2 to 5 do
>   for jj from kk to 5 do
>     R[jj,kk]:= (4^(kk-1)*R[jj,kk-1]-R[jj-1,kk-1])/(4^(kk-1)-1) od od;

> seq([ Je - R[nn,nn], R[nn,nn-1]-R[nn,nn]], nn=2..5);

[-.014159838358561197919, .038187742233276367188],

[-.000346307953198750807, .000863345650335152949],

-5                      -5
[-.1179364820321403*10  , .5392634193412958*10  ],

-9                   -8
[-.388051072754*10  , .4605378004879*10  ]
# In each of these cases the actual error in R[n,n] is less in absolute value than R[n,n-1] - R[n,n].
> f:= x -> x^41: J:= int(f(x), x=a..b); Je:= evalf(J);

J := 1/42

Je := .023809523809523809524
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:
> for kk from 2 to 5 do
>   for jj from kk to 5 do
>     R[jj,kk]:= (4^(kk-1)*R[jj,kk-1]-R[jj-1,kk-1])/(4^(kk-1)-1) od od;

> seq([ Je - R[nn,nn], R[nn,nn-1]-R[nn,nn]], nn=2..5);

[-.059526323670175724408, .041666038130131959938],

[-.015824956803071526056, .002731335429194012396],

[-.001879668742061791717, .000217895125953277100],

-5
[-.000062025551328686623, .7100168713801191*10  ]
# In each of these cases the actual error in R[n,n] is greater in absolute value than R[n,n-1] - R[n,n].
# (c) No, it isn't, as we see above for k = 41.
# (d) R[n,n] is exactly correct for integrating powers up to x^(2*n-1).  So for x^41 we would need
# R[21,21].   This would require calculating the Trapezoid Rule with up
# to 2^21 intervals.
> 2^21;

2097152
# That's not very practical.

```