<jupyter_start><jupyter_text>Use Trapezoid, Simpson's Rule, and Romberg integration<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>Define the function to integrate<jupyter_code>def func(x):
    a = -2.0
    b = 10.0
    
    return (np.exp(a*x))*np.cos(b*x)<jupyter_output><empty_output><jupyter_text>Define the integral to check answer<jupyter_code>def func_integral(x):
    return np.exp(-2*x)*(5*np.sin(10*x)-np.cos(10*x))/52<jupyter_output><empty_output><jupyter_text>Define the core of the trapezoid method<jupyter_code>def trapezoid_core(f, x, h):
    return 0.5*h*(f(x+h) + f(x))<jupyter_output><empty_output><jupyter_text>Define wrapper function to perform trapezoid method<jupyter_code>def trapezoid_method(f,a,b,N):
    
    # f is function to integrate
    # a is lower limit
    # b is upper limit
    # N is the number of chunks to use
    
    # define x values for integral
    x = np.linspace(a,b,N)
    h = x[1]-x[0]
    
    # define the value of the integral
    Fint = 0.0
    
    # perform the integral using the trapezoid method
    for i in range(0,N-1,1):
        Fint += trapezoid_core(f,x[i],h)
    
    # return the answer
    return Fint<jupyter_output><empty_output><jupyter_text>Define the core of Simpson's method<jupyter_code>def simpsons_core(f, x, h):
    return h*( f(x) + 4.*f(x+h) + f(x+2*h))/3.<jupyter_output><empty_output><jupyter_text>Define a wrapper function for Simpson's method<jupyter_code>def simpsons_method(f,a,b,N):
    
    # f is function to integrate
    # a is lower limit
    # b is upper limit
    # N is the number of function evaluations to use
    # note: number of chunks will be N-1
    
    # define x values for integral
    x = np.linspace(a,b,N)
    h = x[1]-x[0]
    
    # define the value of the integral
    Fint = 0.0
    
    # perform the integral using Simpson's method
    for i in range(0,N-2,2):
        Fint += simpsons_core(f,x[i],h)
        
    # apply simpson's over the last interval if N is even
    if ((N%2)==0):
        Fint += simpsons_core(f, x[-2],0.5*h)
    
    # return the answer
    return Fint<jupyter_output><empty_output><jupyter_text>Define the romberg integration core<jupyter_code>def romberg_core(f,a,b,i):
    
    # find the width of the interval
    h = b-a
    
    # set the increment between new function evaluations
    dh = h/2.**(i)
    
    # set the cofactor
    K = h/2.**(i+1)
    
    # sum the function evaluations
    M = 0.0
    for j in range(2**i):
        M += f(a + 0.5*dh + j*dh)
    
    return K*M<jupyter_output><empty_output><jupyter_text>Define a wrapper function for Romberg<jupyter_code>def romberg_integration(f,a,b,tol):
    
    i = 0
    
    # define max number of iterations
    imax = 1000
    
    # Define an error estimate, set to large value
    delta = 100.0*np.fabs(tol)
    
    # set an array of integral answers
    I = np.zeros(imax,dtype=float)
    
    # get zeroth romberg iteration
    I[0] = 0.5*(b-a)* (f(a) + f(b))
    
    # iterate by 1
    i += 1
    
    while(delta>tol):
        
        # find romberg iteration
        I[i] = 0.5*I[i-1] + romberg_core(f,a,b,i)
        
        # compute new fractional error estimate
        delta = np.fabs((I[i] - I[i-1])/I[i])
        
        print (i, I[i],I[i-1],delta)
        if (i >= 20):
            print("This could take a few mintues, please wait..")
        
        if (delta>tol):
            i+=1
            if(i>imax):
                print("max iterations reached")
                raise StopIteration('Stopping iterations after ',i)
    
    return I[i], i<jupyter_output><empty_output><jupyter_text>Perform the integration using each method<jupyter_code>lower_bnd = 0.0
upper_bnd = np.pi
tolerance = 1.0e-6

Answer = func_integral(upper_bnd) - func_integral(lower_bnd)
print("True Answer: ", Answer)



# perfrom Trapezoid integration

# set variable for trapezoid answer to a high value
trap_ans = 10000.

trap_intervals = 1

while (np.fabs(Answer - trap_ans) > tolerance):
    trap_intervals += 1
    trap_ans = trapezoid_method(func, lower_bnd, upper_bnd, trap_intervals)

print("Answer from Trapezoid method: ", trap_ans, " Required intervals: ", trap_intervals)




# perform Simpson's method integration

# set variable for Simpson's method answer to a high value
simpsons_ans = 10000.

simpsons_intervals = 1

while (np.fabs(Answer - simpsons_ans) > tolerance):
    simpsons_intervals += 1
    simpsons_ans = simpsons_method(func, lower_bnd, upper_bnd, simpsons_intervals)

print("Answer from Simpson's method: ", simpsons_ans, " Required intervals: ", simpsons_intervals)



# perform the Romberg method integration

romberg_ans, romberg_intervals = romberg_integration(func, lower_bnd, upper_bnd, tolerance)

print("Answer from Romberg method: ", romberg_ans, " Iterations required: ", romberg_intervals)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11131349929164137, 'fraction_numerical': 0.017607771706132362, 'mean_word_length': 3.1459731543624163, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '124143', 'license': 'MIT', 'n_tokens_mistral': 1876, 'n_tokens_neox': 1772, 'n_words': 448, 'path': 'hw-5.ipynb', 'repo_name': 'natewachter/astr-119-hw-5', 'size': '4941'}