# Some Python code for use at https://sagecell.sagemath.org/ ********** # Numerical Limits def f(x): return (x**2-1)/(x-1) c = 1.0 for i in range(5): e = (10.0)**(-i-1) x1 = c + e print( "x = ", x1, "f(x) = ", f(x1) ) x2 = c - e print( "x = ", x2, "f(x) = ", f(x2) ) print() ********** # Newton's Method def f(x): return x - cos(x) def fp(x): return 1.0 + sin(x) x = 1.0 print(x) for i in range(10): x = x - f(x)/fp(x) print(x) ********** # Simple bisection def f(x): return x - cos(x) a = 0.0 # f(a) and f(b) must be nonzero and have different signs b = 1.0 for i in range(20): c = ( a + b ) / 2 print(c) x = f(c) * f(a) if x == 0: print( "Found it." ) break elif x < 0: b = c else: a = c ********** # Trapzoid Rule def f(x): return 1.0/x a = 1.0 b = 2.0 N = 8 h = ( b - a ) / (1.0*N) x = float(a) s = ( f(a) + f(b) ) / 2.0 for i in range(N-1): x = x + h s = s + f(x) print( "Trapezoid rule approximation:", h * s) ********** # Simpson's rule def f(x): return x**2 + x a = 0.0 b = 3.0 N = 30 # Must be even h = ( b - a ) / (1.0*N) h2 = 2 * h s0 = f(a) + f(b) x = a - h s1 = 0.0 for i in range(0,N,2): x = x + h2 s1 = s1 + f(x) s1 = 4 * s1 x = float(a) s2 = 0.0 for i in range(1,N-1,2): x = x + h2 s2 = s2 + f(x) s2 = 2 * s2 s = h * ( s0 + s1 + s2 ) / 3.0 print( "Simpson's rule approximation:", s )