]> 26.3 Generalization of the Trapezoid Rule

26.3 Generalization of the Trapezoid Rule

Use of the trapezoid rule, which is substantially better than use of the left hand rule for approximating integrals numerically, can be applied here if you can find a way to calculate f ( x , y ) at the right ends of the intervals when you only have an estimate for y at the left end.

An obvious first approximation to doing this is to approximate y at the right end of the interval using the linear approximation to y defined at the left end. The resulting rule is

y ( x + d ) y ( x ) = ( f ( x , y ( x ) ) + f ( x + d , y ( x ) + d * f ( x , y ( x ) ) ) d 2

This rule consists of approximating the difference between the values of y at the ends of the interval by half of d multiplied the sum of the derivative of f at the left end and the linear approximation to the derivative at the right end defined at the left end. When f does not depend on y we get the usual trapezoid rule.

Another way to look at this for subinterval from x to x + d is by defining "iterations" of the left hand rule, as follows

f 0 ( x ) = f ( x , y ( x ) ) f 1 ( x + d ) = f ( x + d , y ( x ) + f 0 ( x ) * d ) f j ( x + d ) = f ( x + d , y ( x ) + f j 1 ( x + d ) * d )

In these terms the computation rule here is

y ( x + d ) = y ( x ) + f 0 ( x ) + f 1 ( x + d ) 2

The left hand rule is off because y changes over the interval. The linear approximation to y applied here is off only because y 's derivative changes, and this is a second derivative effect.

Thus the error it causes is quadratic in the interval size, and is on a par with the error intrinsic to the trapezoid rule. We therefore expect this rule to give results that improve in accuracy by a factor of four when N is doubled.

Again, there is no great complication in setting up a spreadsheet to compute the predictions of this rule for any N and you can extrapolate as before with it.

It has the advantage that you can start with the factor 4 extrapolation because accuracy improvement by a factor of 4 on doubling the number of points is built into its structure.

The rule no longer has the same weight structure as the trapezoid rule because when you compute f ( x , y ) at a given intermediate point from the left you are using the linear approximation at the previous point, while when you compute it in the interval to its right you use the value of y computed from the rule itself in the previous interval. Such is life.

Here are the entries in rows 9 and 10 that can be used in place of those above to produce this computation, when those in red are copied down. The red entries in columns D and E are given for the differential equation y ' = x + y . They have to be changed and the results copied down, in order to switch to a different differential equation.

B

C

D

E

 

X

y = y ( x d ) + ( f 0 ( x d ) + f 1 ( x ) ) * d 2 f 0 ( x ) = x + y ( x ) f 1 ( x + d ) = x + d + y ( x ) + d * f 0 ( x )

8

=B3

=B6

=B9+C9

=B10+C9+(B10-B9)*D9

9

=B9+$B$7

=C9+(D9+E9)*(B10-B9)/2

=B10+C10

=B12+C10+(B12-B10)*D10

10

Exercise 26.2 Compare the results obtained using the rule above, with those obtained using the left hand rule, upon the same problem.

Here are results obtained using this trapezoid like method with various levels of extrapolation

N\L#

L 1 L 2 L 3 L 4 L 5 L 6 L 7 L 8

1

5

             

2

9.5

11

           

4

10.9458

11.42773

11.48883929

         

8

11.52449

11.71739

11.75877194

11.77676745

       

16

11.70817

11.76939

11.77681781

11.77802087

11.7780613

     

32

11.75976

11.77696

11.77804039

11.77812189

11.77812515

11.77812617

   

64

11.77341

11.77796

11.77810835

11.77811289

11.77811259

11.7781124

11.77811229

 

128

11.77692

11.77809

11.77811199

11.77811223

11.77811221

11.7781122

11.7781122

11.7781122

The accuracy of these calculations for this problem are shown in the following table

N\L#

1

2

3

4

5

6

7

8

1

-0.57548

 

           

2

-0.19342

-0.06606

           

4

-0.07067

-0.02975

-0.02456

 

 

     

8

-0.02153

-0.00516

-0.00164

-0.00011

 

     

16

-0.00594

-0.00074

-0.00011

-7.8E-06

-4E-06

     

32

-0.00156

-9.8E-05

-6.1E-06

8.23E-07

1.1E-06

1.19E-06

 

 

64

-0.0004

-1.3E-05

-3.3E-07

5.84E-08

3.4E-08

1.68E-08

7.6E-09

 

128

-0.0001

-1.6E-06

-1.8E-08

2.51E-09

7.1E-10

1.84E-10

5.3E-11

2.4E-11

Proportional Error for f 0 ( x ) + f 1 ( x + d ) 2 Rule y ' = y + z , y ( 0 ) = 1 , finding y ( 2 )

   

You will note that the estimates here without extrapolation are a little better than those for the first iteration of the previous method, by a factor of 6 for N = 128 . However, after extrapolation the results are more accurate by a factor of thousands than in the previous table. See First Order ODE Applet