First, we study richardson extrapoloation by reminding ourselves how the trapezoidal rule worked: >> for i=1:9, [area,h,error(i)] = trap(1,2,2^i+1); end >> for i=1:8, ratio(i) = error(i)/error(i+1); end >> error error = Columns 1 through 6 3.0256e-003 7.6109e-004 1.9058e-004 4.7665e-005 1.1917e-005 2.9794e-006 Columns 7 through 9 7.4486e-007 1.8622e-007 4.6554e-008 >> ratio The ratio's going to 4. ratio = Columns 1 through 6 3.9753e+000 3.9935e+000 3.9984e+000 3.9996e+000 3.9999e+000 4.0000e+000 Columns 7 through 8 4.0000e+000 4.0000e+000 Now we try the Richardson extrapolation: >> for i=1:9, [area,h,error(i)] = richardson_trap(1,2,2^i+1); end >> for i=1:8, ratio(i) = error(i)/error(i+1); end >> error error = Columns 1 through 6 6.2596e-006 4.1111e-007 2.6047e-008 1.6337e-009 1.0219e-010 6.3880e-012 Columns 7 through 9 3.9946e-013 2.4869e-014 4.2188e-015 >> ratio The ratio's going to 16. ratio = Columns 1 through 6 1.5226e+001 1.5783e+001 1.5944e+001 1.5986e+001 1.5998e+001 1.5992e+001 Columns 7 through 8 1.6063e+001 5.8947e+000 Now we do Simpsons rule: >> for i=1:9, [area,h,error(i)] = simp(1,2,2^i+1); end >> for i=1:8, ratio(i) = error(i)/error(i+1); end >> error error = Columns 1 through 6 8.5909e-005 6.2596e-006 4.1111e-007 2.6047e-008 1.6337e-009 1.0219e-010 Columns 7 through 9 6.3884e-012 4.0012e-013 2.5979e-014 >> ratio The ratio's going to 16. ratio = Columns 1 through 6 1.3724e+001 1.5226e+001 1.5783e+001 1.5944e+001 1.5986e+001 1.5997e+001 Columns 7 through 8 1.5966e+001 1.5402e+001 Now we try Richardson extrapolation on Simpsons rule: >> for i=1:9, [area,h,error(i)] = richardson_simp(1,2,2^i+1); end >> for i=1:8, ratio(i) = error(i)/error(i+1); end >> error error = Columns 1 through 6 9.4972e-007 2.1208e-008 3.7638e-010 6.1020e-012 9.5479e-014 1.5543e-015 Columns 7 through 9 1.1102e-015 1.1102e-015 -4.4409e-016 >> ratio The ratio appears to be going to 64 before it hits round-off error ratio = Columns 1 through 6 4.4782e+001 5.6347e+001 6.1681e+001 6.3909e+001 6.1429e+001 1.4000e+000 Columns 7 through 8 1.0000e+000 -2.5000e+000 Now we turn to numerical differentiation: >> for i=1:9, err(i) = 1/2*2^(-1/2) - d_right(2,1/2^i); end >> err err = Columns 1 through 6 1.9703e-002 1.0408e-002 5.3581e-003 2.7198e-003 1.3704e-003 6.8785e-004 Columns 7 through 9 3.4459e-004 1.7247e-004 8.6275e-005 >> for i=1:8, ratio(i) = err(i)/err(i+1); end >> ratio The ratio's going to 2, as expected. ratio = Columns 1 through 6 1.8931e+000 1.9424e+000 1.9700e+000 1.9847e+000 1.9923e+000 1.9961e+000 Columns 7 through 8 1.9981e+000 1.9990e+000 Now we try the midpoint rule >> for i=1:9, err(i) = 1/2*2^(-1/2) - d_mid(2,1/2^i); end >> err err = Columns 1 through 6 -2.8406e-003 -6.9530e-004 -1.7293e-004 -4.3177e-005 -1.0791e-005 -2.6975e-006 Columns 7 through 9 -6.7435e-007 -1.6859e-007 -4.2147e-008 >> for i=1:8, ratio(i) = err(i)/err(i+1); end >> ratio The ratio's going to 4, as expected. ratio = Columns 1 through 6 4.0854e+000 4.0207e+000 4.0051e+000 4.0013e+000 4.0003e+000 4.0001e+000 Columns 7 through 8 4.0000e+000 4.0000e+000 Now we try a program where we use the trapezoidal rule to integrate a function that involves a derivative. The derivative is approximated using the right-hand rule. >> for i=1:9, goof(1,2,2^i+1); end >> for i=1:9, area(i) = goof(1,2,2^i+1); end >> area area = Columns 1 through 6 3.9385e-002 3.8346e-002 3.8624e-002 3.9014e-002 3.9281e-002 3.9434e-002 Columns 7 through 9 3.9516e-002 3.9557e-002 3.9579e-002 >> for i=1:7, ratio(i) = (area(i)-area(i+1))/(area(i+1)-area(i+2)); end >> ratio The ratio's going to 2, as expected ratio = Columns 1 through 6 -3.7452e+000 7.0962e-001 1.4630e+000 1.7488e+000 1.8779e+000 1.9397e+000 Column 7 1.9701e+000 Now we repeat the exercise, but we approximate the derivative using the midpoint rule. >> for i=1:9, area(i) = goof(1,2,2^i+1); end >> for i=1:7, ratio(i) = (area(i)-area(i+1))/(area(i+1)-area(i+2)); end >> ratio Now the ratio is going to 4, as desired. ratio = Columns 1 through 6 4.4659e+000 4.1050e+000 4.0256e+000 4.0064e+000 4.0016e+000 4.0004e+000 Column 7 4.0001e+000 diary off