I lost the diary from class 1/13/99 when I tripped over the power-cord. But here's a variant of it. >> diary 1_13_99 >> help trap this has the integration for the trapezoidal rule, n can be an even or odd number. function [area,h,error] = trap(a,b,n) I've chosen f(x)=cos(x) and f_int(x)=sin(x). I'm integrating from 1 to 3 using the trapezoidal rule. As discussed in class, we expect the errors to decrease like h^2. >> for i=1:8, [area(i),h(i),error(i)] = trap(1,3,3^i+1); end >> format short e To five sig figs, the areas have converged. If I'd looked at the area after typing "format long e" I would see more sig figs and see that the area hasn't converged to 16 sig figs. >> area area = Columns 1 through 6 -6.7422e-01 -6.9747e-01 -7.0003e-01 -7.0032e-01 -7.0035e-01 -7.0035e-01 Columns 7 through 8 -7.0035e-01 -7.0035e-01 The error is has not reached round-off level yet. (for matlab on hans that's around 10^(-15). >> error error = Columns 1 through 6 -2.6133e-02 -2.8845e-03 -3.2026e-04 -3.5582e-05 -3.9535e-06 -4.3928e-07 Columns 7 through 8 -4.8809e-08 -5.4232e-09 To see if the errors are decreasing like h^2, we can construct two ratios. One from the errors and one from the computed areas. (We often don't know the exact errors computationally, but we always know the computed areas.) >> for i=1:7, ratio1(i)=error(i)/error(i+1); end >> for i=1:6, ratio2(i)=(area(i)-area(i+1))/(area(i+1)-area(i+2)); end Since I was dividing h by 3 every time, the ratios should go to 9. >> ratio1 ratio1 = Columns 1 through 6 9.0599e+00 9.0066e+00 9.0007e+00 9.0001e+00 9.0000e+00 9.0000e+00 Column 7 9.0000e+00 >> ratio2 ratio2 = 9.0666e+00 9.0073e+00 9.0008e+00 9.0001e+00 9.0000e+00 9.0000e+00 Recover the rate of convergence by taking the logarithms. We see the power 2 as desired. >> log(ratio1)/log(3) ans = Columns 1 through 6 2.0060e+00 2.0007e+00 2.0001e+00 2.0000e+00 2.0000e+00 2.0000e+00 Column 7 2.0000e+00 >> log(ratio2)/log(3) ans = 2.0067e+00 2.0007e+00 2.0001e+00 2.0000e+00 2.0000e+00 2.0000e+00 What happens when we run into round-off error? I test this by integrating from 1 to 1.01 instead of 1 to 3. (Since I don't want to get slowed down by trying to do the above computation but for 3^15 intervals.) >> for i=1:8, [area(i),h(i),error(i)] = trap(1,1.01,3^i+1); end >> for i=1:7, ratio1(i)=error(i)/error(i+1); end >> for i=1:6, ratio2(i)=(area(i)-area(i+1))/(area(i+1)-area(i+2)); end The error's hit round-off: >> error error = Columns 1 through 6 4.9638e-09 5.5153e-10 6.1281e-11 6.8090e-12 7.5655e-13 8.4060e-14 Columns 7 through 8 9.3398e-15 1.0400e-15 The ratios go to 9 and then move away when we get close to round-off. >> format long e >> ratio1 ratio1 = Columns 1 through 3 9.000001486154204e+00 8.999999731076866e+00 9.000003312001917e+00 Columns 4 through 6 9.000005732314662e+00 9.000134138162307e+00 9.000278603268946e+00 Column 7 8.980817347789825e+00 >> ratio1-9 ans = Columns 1 through 3 1.486154204144441e-06 -2.689231344277232e-07 3.312001917166185e-06 Columns 4 through 6 5.732314662054705e-06 1.341381623074511e-04 2.786032689456874e-04 Column 7 -1.918265221017457e-02 >> ratio2 ratio2 = Columns 1 through 3 9.000001705538878e+00 8.999999283461420e+00 9.000003009463041e+00 Columns 4 through 6 8.999989681852830e+00 9.000116080652838e+00 9.002717107325740e+00 >> ratio2-9 ans = Columns 1 through 3 1.705538878127300e-06 -7.165385795815382e-07 3.009463041436788e-06 Columns 4 through 6 -1.031814717045165e-05 1.160806528375247e-04 2.717107325739931e-03 >> format short e Now we repeat the exercise with Simpsons rule. >> help simp this has the integration for Simpsons rule, n must be odd function [area,h,error] = simp(a,b,n) >> for i=1:12, [area(i),h(i),error(i)] = simp(1,3,2^i+1); end >> for i=1:10, ratio2(i)=(area(i)-area(i+1))/(area(i+1)-area(i+2)); end >> for i=1:11, ratio1(i)=error(i)/error(i+1); end >> error error = Columns 1 through 6 4.4082e-03 2.5061e-04 1.5312e-05 9.5168e-07 5.9397e-08 3.7110e-09 Columns 7 through 12 2.3192e-10 1.4495e-11 9.0605e-13 5.7510e-14 2.2204e-15 5.5511e-16 The ratios should go to 16. Which they do. But they move away when we get close to round-off. >> ratio1 ratio1 = Columns 1 through 6 1.7590e+01 1.6366e+01 1.6090e+01 1.6022e+01 1.6006e+01 1.6001e+01 Columns 7 through 11 1.6000e+01 1.5998e+01 1.5755e+01 2.5900e+01 4.0000e+00 >> ratio2 ratio2 = Columns 1 through 6 1.7670e+01 1.6385e+01 1.6094e+01 1.6023e+01 1.6006e+01 1.6002e+01 Columns 7 through 10 1.6000e+01 1.6015e+01 1.5347e+01 3.3200e+01 The logarithms should go to 4... >> log(ratio1)/log(2) ans = Columns 1 through 6 4.1367e+00 4.0327e+00 4.0081e+00 4.0020e+00 4.0005e+00 4.0001e+00 Columns 7 through 11 4.0000e+00 3.9998e+00 3.9777e+00 4.6949e+00 2.0000e+00 >> log(ratio2)/log(2) ans = Columns 1 through 6 4.1432e+00 4.0343e+00 4.0085e+00 4.0021e+00 4.0005e+00 4.0001e+00 Columns 7 through 10 4.0000e+00 4.0013e+00 3.9399e+00 5.0531e+00 Now I do a case where I've artificially introduced a bug. Specifically, I added "area = area + .1*f(x(n))" before the line "area = area*h/3" >> for i=1:12, [area(i),h(i),error(i)] = simp(1,3,2^i+1); end >> for i=1:11, ratio1(i)=error(i)/error(i+1); end >> for i=1:10, ratio2(i)=(area(i)-area(i+1))/(area(i+1)-area(i+2)); end The error is still decreasing: the patient lives, but the cleverness of simpsons rules is wasted. The error is decreasing no more quickly than it would for a random Riemann sum. In fact, the error is like C h^2 + D h. If the bug had a smaller coefficient, if I'd added .001*f(x(n)) then the effect would be that D is smaller. This means that the error would look like h^2 until h got much smaller and I might fool myself. You can see a ghost of this in the ratios. The first ones are close to 4. >> error error = Columns 1 through 6 7.7082e-03 1.9006e-03 8.4031e-04 4.1345e-04 2.0631e-04 1.0313e-04 Columns 7 through 12 5.1562e-05 2.5781e-05 1.2891e-05 6.4453e-06 3.2226e-06 1.6113e-06 >> ratio1 ratio1 = Columns 1 through 6 4.0557e+00 2.2618e+00 2.0324e+00 2.0040e+00 2.0005e+00 2.0001e+00 Columns 7 through 11 2.0000e+00 2.0000e+00 2.0000e+00 2.0000e+00 2.0000e+00 >> ratio2 ratio2 = Columns 1 through 6 5.4773e+00 2.4839e+00 2.0607e+00 2.0076e+00 2.0009e+00 2.0001e+00 Columns 7 through 10 2.0000e+00 2.0000e+00 2.0000e+00 2.0000e+00 Now I've removed the bug and will try some other functions. First I do f(x)=x^(1/2), integraing from 0 to 1. Since f'''' is infinite at 0, I don't have any right for the errors to decrease like h^4. >> for i=1:12, [area(i),h(i),error(i)] = simp(0,1,2^i+1); end >> for i=1:11, ratio1(i)=error(i)/error(i+1); end The errors are decreasing... >> error error = Columns 1 through 6 2.8595e-02 1.0140e-02 3.5874e-03 1.2685e-03 4.4848e-04 1.5856e-04 Columns 7 through 12 5.6061e-05 1.9820e-05 7.0076e-06 2.4776e-06 8.7595e-07 3.0969e-07 They're decreasing like h^(3/2): >> log(ratio1)/log(2) ans = Columns 1 through 6 1.4957e+00 1.4991e+00 1.4998e+00 1.5000e+00 1.5000e+00 1.5000e+00 Columns 7 through 11 1.5000e+00 1.5000e+00 1.5000e+00 1.5000e+00 1.5000e+00 Now I've taken f(x)=x^(3/2) >> for i=1:12, [area(i),h(i),error(i)] = simp(0,1,2^i+1); end >> for i=1:11, ratio1(i)=error(i)/error(i+1); end The errors are decreasing like h^(5/2): >> log(ratio1)/log(2) ans = Columns 1 through 6 2.4554e+00 2.4832e+00 2.4939e+00 2.4978e+00 2.4992e+00 2.4997e+00 Columns 7 through 11 2.4999e+00 2.5000e+00 2.5000e+00 2.5000e+00 2.5001e+00 Now I've taken f(x)=x^(5/2) >> for i=1:12, [area(i),h(i),error(i)] = simp(0,1,2^i+1); end >> for i=1:11, ratio1(i)=error(i)/error(i+1); end The errors are decreasing like h^(7/2). Although it's less convincing here. It'd be more convincing if we could do the computation in quadruple precision. >> log(ratio1)/log(2) ans = Columns 1 through 6 3.2969e+00 3.3672e+00 3.4120e+00 3.4407e+00 3.4595e+00 3.4720e+00 Columns 7 through 11 3.4805e+00 3.4864e+00 3.4905e+00 3.5016e+00 3.3158e+00 >> error error = Columns 1 through 6 1.1965e-03 1.2174e-04 1.1798e-05 1.1084e-06 1.0208e-07 9.2799e-09 Columns 7 through 12 8.3631e-10 7.4925e-11 6.6853e-12 5.9480e-13 5.2514e-14 5.2736e-15 >> diary off