# function difference_table(f,fprime,x) # Accompanying program for the text # # Classical and Modern Numerical Analysis: # Theory, Methods and Practice # by Azmy S. Ackleh, Edward J. Allen, # R. Baker Kearfott, and Padmanabhan Seshaiyer # # (Taylor and Francis / CRC Press, 2009) # # Issuing the command # difference_table('f','fprime',x) # causes a table of difference quotients to be formed using the # difference_quotient function. The string f can be either the name of a # Matlab intrinsic function or a user-supplied m-file, while fprime is a # similar pointer to its derivative. The step sizes used are h = 1/4, # h = 1/8, ... h = 1/4^{30}. For example, issuing the command # # difference_table('sin','cos',pi/4) # # produces the followiong table: # # 1 2.50e-001 0.6118351194488110 -9.53e-002 # 2 6.25e-002 0.6845566203276636 -2.26e-002 # 3 1.56e-002 0.7015538499518499 -5.55e-003 # 4 3.91e-003 0.7057239167465070 -1.38e-003 # 5 9.77e-004 0.7067614018394579 -3.45e-004 # 6 2.44e-004 0.7070204574170020 -8.63e-005 # 7 6.10e-005 0.7070852015622222 -2.16e-005 # 8 1.53e-005 0.7071013863678672 -5.39e-006 # 9 3.81e-006 0.7071054324915167 -1.35e-006 # 10 9.54e-007 0.7071064440533519 -3.37e-007 # 11 2.38e-007 0.7071066969074309 -8.43e-008 # 12 5.96e-008 0.7071067616343498 -1.96e-008 # 13 1.49e-008 0.7071067765355110 -4.65e-009 # 14 3.73e-009 0.7071067988872528 1.77e-008 # 15 9.31e-010 0.7071068286895752 4.75e-008 # 16 2.33e-010 0.7071070671081543 2.86e-007 # 17 5.82e-011 0.7071075439453125 7.63e-007 # 18 1.46e-011 0.7071075439453125 7.63e-007 # 19 3.64e-012 0.7071228027343750 1.60e-005 # 20 9.09e-013 0.7071533203125000 4.65e-005 # 21 2.27e-013 0.7075195312500000 4.13e-004 # 22 5.68e-014 0.7070312500000000 -7.55e-005 # 23 1.42e-014 0.7109375000000000 3.83e-003 # 24 3.55e-015 0.7187500000000000 1.16e-002 # 25 8.88e-016 0.7500000000000000 4.29e-002 # 26 2.22e-016 1.0000000000000000 2.93e-001 # 27 5.55e-017 0.0000000000000000 -7.07e-001 # 28 1.39e-017 0.0000000000000000 -7.07e-001 # 29 3.47e-018 0.0000000000000000 -7.07e-001 # 30 8.67e-019 0.0000000000000000 -7.07e-001 # # You could also issue, say, the command # # difference_table('xsqm2', 'xsqm2_prime', 1) # # where xsqm2 and xsqm2_prime are supplied by you. (See the in-line # documentation for bisect.m.) from difference_quotient import difference_quotient def difference_table(f,fprime,x): for i in range(1,31): h=4**(-i) value = difference_quotient(f,x,h) error = value - fprime(x) print('{:3d} {:12.2e} {:20.16f} {:12.2e}'.format(i, h, value, error)) if __name__ == '__main__': # test it from numpy import sin, cos, pi difference_table( sin, cos, pi/4 ) ''' $ py difference_table.py 1 2.50e-01 0.6118351194488110 -9.53e-02 2 6.25e-02 0.6845566203276636 -2.26e-02 3 1.56e-02 0.7015538499518499 -5.55e-03 4 3.91e-03 0.7057239167465070 -1.38e-03 5 9.77e-04 0.7067614018394579 -3.45e-04 6 2.44e-04 0.7070204574170020 -8.63e-05 7 6.10e-05 0.7070852015622222 -2.16e-05 8 1.53e-05 0.7071013863678672 -5.39e-06 9 3.81e-06 0.7071054324915167 -1.35e-06 10 9.54e-07 0.7071064440533519 -3.37e-07 11 2.38e-07 0.7071066969074309 -8.43e-08 12 5.96e-08 0.7071067616343498 -1.96e-08 13 1.49e-08 0.7071067765355110 -4.65e-09 14 3.73e-09 0.7071067988872528 1.77e-08 15 9.31e-10 0.7071068286895752 4.75e-08 16 2.33e-10 0.7071070671081543 2.86e-07 17 5.82e-11 0.7071075439453125 7.63e-07 18 1.46e-11 0.7071075439453125 7.63e-07 19 3.64e-12 0.7071228027343750 1.60e-05 20 9.09e-13 0.7071533203125000 4.65e-05 21 2.27e-13 0.7075195312500000 4.13e-04 22 5.68e-14 0.7070312500000000 -7.55e-05 23 1.42e-14 0.7109375000000000 3.83e-03 24 3.55e-15 0.7187500000000000 1.16e-02 25 8.88e-16 0.7500000000000000 4.29e-02 26 2.22e-16 1.0000000000000000 2.93e-01 27 5.55e-17 0.0000000000000000 -7.07e-01 28 1.39e-17 0.0000000000000000 -7.07e-01 29 3.47e-18 0.0000000000000000 -7.07e-01 30 8.67e-19 0.0000000000000000 -7.07e-01 '''