<jupyter_start><jupyter_text>HW 5: Numerical Integration<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def func(x):
    return np.exp(-2*x) * np.cos(10*x)<jupyter_output><empty_output><jupyter_text>Trapazoidal Method<jupyter_code>def trapazoid(lower_bound, upper_bound, intervals):                #Function performing Trapazoid method
    x1 = lower_bound
    x2 = upper_bound
    n = intervals

    delta_x = (x2 - x1) / n                                        #Trapazoid width
    i = 0                                                          #Iteration tracker
    total = 0.0                                                    #Running total of integral
    
    #print(delta_x)
    
    for i in range (n):
        low = x1 + (i*delta_x)
        high = x1 + ((i+1)*delta_x)
        #print((delta_x/2) * (func(low) + func(high)))
        total += (delta_x/2) * (func(low) + func(high))
        #print(total)
        #print(low)
        #print(high)
        
    return total<jupyter_output><empty_output><jupyter_text>Simpson's<jupyter_code>def simpson(lower_bound, upper_bound, intervals):
    x1 = lower_bound
    x2 = upper_bound
    n = intervals
    
    if(n%2 == 1):
        print("Odd Interval")
        return 0
    
    delta_x = (x2 - x1) / n
    i = 0
    total = 0.0                              #Running total
    
    #print(delta_x)
    
    for i in range (n):
        if(i == 0):                          #First term
            total += func(x1)
            #print(total)
            
        elif(i == n-1):
            total += func(x2)
            #print("Loop finished")
            #print(total)
            
        elif(i%2 == 0):                      #Even Term *4
            total += 4*func(x1+i*delta_x)
            #print(i)
            #print("^*4")
            #print(total)
            
        elif(i%2 == 1):                      #Odd Term *2
            total += 2*func(x1+i*delta_x)
            #print(i)
            #print(total)
            
    total *= delta_x / 3
    
    return total<jupyter_output><empty_output><jupyter_text>Romberg<jupyter_code>def base(a, h, i):
    dh = h/2.**(i)
    running = 0.0
    
    for j in range(2**i):
        running += func(a + 0.5*dh + j*dh)
    
    running *= (dh*0.5)
    
    return running
    

def romberg(lower_bound, upper_bound):
    a = lower_bound
    b = upper_bound
    h = b-a
    imax = 30
    tol = 1.0e-6
    I = np.zeros(imax, dtype=float)
    d = 100.
    
    for i in range (0, imax):
        if(i == 0):                                  #Set zeroth Romberg
            #print("Test")
            I[i] = 0.5 * h * (func(a) + func(b))
        
        else:
            I[i] = (0.5 * I[i-1]) + base(a, h, i)
            d = np.fabs((I[i]-I[i-1])/I[i])
            #print("d is: ", d)
            
            if(d < tol):
                return I[i]
    
    print("Max iterations reached, integration is a best estimate.")
    return I[imax]
ans = 0.019195                               #Used my calculator wich does definite integrals for answer
tol = 1.0e-6

for i in range (1, 2000):
    if(np.fabs(ans - trapazoid(0, np.pi, i)) < tol):
        print("Using Trapazoid Method integral is: ", trapazoid(0, np.pi, i), "found after: ", i, " iterations." )
        break

for i in range (1, 10000):
    if(np.fabs(ans - simpson(0, np.pi, 2*i)) < 1.0e-3):
        print("Using Simpson's Method integral is: ", simpson(0, np.pi, i), "found after: ", i, " iterations." )
        break

print("Simpson's used lower accuracy, didn't converge easily with 1e-6.")

print("Using Romberg method: ", romberg(0, np.pi))

print("Finished")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13050526884625777, 'fraction_numerical': 0.02647932991083491, 'mean_word_length': 1.5654885654885655, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '397181', 'license': 'MIT', 'n_tokens_mistral': 1317, 'n_tokens_neox': 1250, 'n_words': 297, 'path': 'hw-5.ipynb', 'repo_name': 'nbkroeker/astr-119-hw-5', 'size': '3701'}