<jupyter_start><jupyter_text>Assignment 8 Problem 1.1We will investigate the following numerical methods for estimating integrals: Simpson, adaptive trapezoid, and traditional trapezoid.Both the trapezoid and Simpson methods divide the x-range over which a function is integrated into subdivisions. To use both method, you have to specify the number of subdivisions you want beforehand. The trapezoid method uses a trapezoid to evaluate the area within each subdivision and then add them up. In contrast, the Simpson method uses a polynomial to match the shape of the curve, and evaluate the area within each subdivision, as a result, we expect the Simpson method to give us more accurate result. The adaptive trapezoid method uses the same approach as the traditional trapezoid method, but instead of specifying the number of subdivisions beforehand, you specify a level of accuracy you want, and the adaptive method will divide the x-range into more and more subdivisions until that level of accuracy is reached. Let's see them in action.First of all, I want to enable linting. I typed the following commands:%load_ext pycodestyle_magic%pycodestyle_oBut there is an error: "pycodestyle module not found".Turns out, pycodestyle_magic module has not been installed yet. So I installed pycodestyle_magic first with: %pip install pycodestyle_magicAnd I ran the command again. There is another error: "flake8 module not found".So, in order to enable linting, I need to install both modules, flake8 and pycodestyle_magic.I did so, and I could enable linting finally as shown below. :)<jupyter_code>%pip install flake8 pycodestyle_magic
# This is to enable linting, making sure we write our python code in style.
%load_ext pycodestyle_magic
%pycodestyle_on
# I wonder if I can import all three methods in a single command?
# Apparently I can.
from integrals import adaptive_trapezoid, trapezoid, simpson

# This magical command imports both numpy and matplotlib in one go.
%pylab inline
# adaptive_trapezoid(the function, lower limit, upper limit, level of accuracy, output=True or false)
# simpson(the function, lower limit, upper limit, the number of subdivisions)
# trapezoid (the function, lower limit, upper limit, the number of subdivisions)
# The above comments are too long, but they specifies the inputs.<jupyter_output>1:80: E501 line too long (101 > 79 characters)
3:80: E501 line too long (80 > 79 characters)<jupyter_text>We first invetigate the function $y = e^x$ in the following range [0, 10].<jupyter_code>calc_e_adaptive = adaptive_trapezoid(np.exp, a=0, b=10, acc=0.5, output=True)<jupyter_output>N = 2,  Integral = 110137.32897403359
N = 2.0,  Integral = 55810.73028252968
N = 4.0,  Integral = 32455.927412306755
N = 8.0,  Integral = 24821.354153296994
N = 16.0,  Integral = 22737.81578168463
N = 32.0,  Integral = 22204.418433198334
N = 64.0,  Integral = 22070.25849666271
N = 128.0,  Integral = 22036.667386590423
N = 256.0,  Integral = 22028.26640638897
N = 512.0,  Integral = 22026.16596105638
N = 1024.0,  Integral = 22025.64083720381<jupyter_text>Here I specified the accuracy to be 0.5, and compared to the actual value of the integral below, indeed it's accurate to within 0.5. Five digits are correct!<jupyter_code>np.exp(10)-np.exp(0)<jupyter_output><empty_output><jupyter_text>As we explained earlier, the adaptive method tries more and more subdivisions until the level of accuracy we specified is reached. In this case it took 11 steps for it to reach that level of accuracy. For the other two methods, trapezoid and Simpson, there is really just 1 step! We specify the number of subdivisons and that's it, the methods then calculate a value for us, done. (It would probably make more sense to compare the efficiency between adaptive trapezoid and adaptive Simpson.)<jupyter_code>calc_e_trapezoid = trapezoid(np.exp, 0, 10, 1024)
print(calc_e_trapezoid)<jupyter_output>22025.640837203788<jupyter_text>As we can see, just one step, and boom we have the answer. Here I specified the number of subdivisions to be 1024. I expected that both the adaptive and traditional trepazoid methods should give me similar results if I make the number of subdivision the same. Since they're essentially the same method. And indeed they are exactly the same answers!<jupyter_code>calc_e_simpson = simpson(np.exp, 0, 10, 1024)
print(calc_e_simpson)<jupyter_output>22025.465795919594<jupyter_text>The Simpson method is very accurate as expected. It uses 1024 subdivisons, and it is accurate up to 6 digits. It is more accurate than the adaptive and traditional trapezoid methods, which produced accurate answers up to 5 digits. I guess the adaptive Simpson method might also be more efficient, since it can more accurately simulate the shape of the curve, so I guess it will take fewer steps for adaptive Simpson to reach the accuracy we want. But we don't have an adaptive Simpson method code to test this hypothesis yet.  The integral we're trying to determine is plotted below:<jupyter_code>a = 0
b = 10
x = np.linspace(a, b, 10000)
y = np.exp(x)
plt.xlabel("x")
plt.ylabel("exp(x)")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()<jupyter_output><empty_output><jupyter_text>---Now we will try different ranges, first x = [0, 2𝜋]. The result we expect should be:<jupyter_code>np.exp(2*np.pi)-np.exp(0)
a = 0
b = 2*np.pi
x = np.linspace(a, b, 10000)
y = np.exp(x)
plt.xlabel("x")
plt.ylabel("exp(x)")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()<jupyter_output><empty_output><jupyter_text>Next, let's do the integrals numerically.<jupyter_code>calc_e_adaptive2 = adaptive_trapezoid(np.exp, a=0, b=2*np.pi, acc=0.5, output=True)
calc_e_trapezoid2 = trapezoid(np.exp, 0, 2*np.pi, 128)
print(calc_e_trapezoid2)
calc_e_simpson2 = simpson(np.exp, 0, 2*np.pi, 128)
print(calc_e_simpson2)<jupyter_output>534.4916727602475<jupyter_text>Both adaptive and traditional trepazoid methods are accurate up to 3 digits, while the Simpson method is accurate up to 7 digits!  ---Next we change the range to x = [-1. 1]. The expected result is:<jupyter_code>np.exp(1)-np.exp(-1)
a = -1
b = 1
x = np.linspace(a, b, 10000)
y = np.exp(x)
plt.xlabel("x")
plt.ylabel("exp(x)")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()
calc_e_adaptive3 = adaptive_trapezoid(np.exp, a=-1, b=1, acc=0.05, output=True)
calc_e_trapezoid3 = trapezoid(np.exp, -1, 1, 8)
print(calc_e_trapezoid3)
calc_e_simpson3 = simpson(np.exp, -1, 1, 8)
print(calc_e_simpson3)<jupyter_output>2.3504530172422795<jupyter_text>The Simpson method is accurate up to 5 digits. While the trapezoid methods are accurate up to 2 digits. Here you might notice that we only used 8 subdivisions to reach the required level of accuracy. This is because the real area is quite small. Still, the Simpson method can give us surprisingly good result with so few subdivisions. --- Problem 1.2We repeat the same steps as above for a new function $y = x^3$. In the range x = [0,10], the result we expect is:<jupyter_code>np.power(10, 4)/4
a = 0
b = 10
x = np.linspace(a, b, 10000)
y = np.power(x, 3)
plt.xlabel("x")
plt.ylabel("x^3")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()
# We define the function first. 
# Somehow we could't put np.power(x, 3) directly into the argument.

def f(x) :
    return np.power(x, 3)
calc_e_adaptive4 = adaptive_trapezoid(f, a=0, b=10, acc=0.5, output=True)
calc_e_trapezoid4 = trapezoid(f, 0, 10, 128)
print(calc_e_trapezoid4)
calc_e_simpson4 = simpson(f, 0, 10, 128)
print(calc_e_simpson4)<jupyter_output>2500.0<jupyter_text>The trepazoid methods are accurate up to 4 digits. The Simpson method gives us the exact answer! Initially it is surprising, but if we think about it, the Simpson method uses a polynomial to match the curve, and in this case, it can find a perfect match since the curve is a polynomial function. ---Changing the range to x = [0, 2𝜋]. The answer we expect is:<jupyter_code>np.power(np.pi*2, 4)/4
a = 0
b = np.pi*2
x = np.linspace(a, b, 10000)
y = np.power(x, 3)
plt.xlabel("x")
plt.ylabel("x^3")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()
calc_e_adaptive5 = adaptive_trapezoid(f, a=0, b=np.pi*2, acc=0.5, output=True)
calc_e_trapezoid5 = trapezoid(f, 0, np.pi*2, 64)
print(calc_e_trapezoid5)
calc_e_simpson5 = simpson(f, 0, np.pi*2, 64)
print(calc_e_simpson5)<jupyter_output>389.63636413600966<jupyter_text>Simpson method gives exact answer. But does it always give us the exact answer when the function is a polynomial? Not necessarily, as We shall see in the next section. I don't know the reason.The trepazoid methods are accurate up to 3 digits. ---Changing the range to x = [-1, 1]. The answer we expect is:<jupyter_code>np.power(1, 4)/4 - np.power(-1, 4)/4
a = -1
b = 1
x = np.linspace(a, b, 10000)
y = np.power(x, 3)
plt.xlabel("x")
plt.ylabel("x^3")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()
calc_e_adaptive6 = adaptive_trapezoid(f, a=-1, b=1, acc=0.05, output=True)
calc_e_trapezoid6 = trapezoid(f, -1, 1, 2)
print(calc_e_trapezoid6)
calc_e_simpson6 = simpson(f, -1, 1, 2)
print(calc_e_simpson6)<jupyter_output>0.0<jupyter_text>All three methods are exact. But it is a coincidence I think. Because the function is odd, when N =2 for the three methods, the function is divided equally across the origin, both areas are equal in magnitude but opposite in sign. This gives us the exact answer 0. For the adaptive trepazoid method, N = 2 is a good enough attempt, it has reached the accuracy level we want. So it stops there. Let's give the traditional trepazoid and Simpson methods more subdivisions, and see what we get:<jupyter_code>calc_e_trapezoid6_2 = trapezoid(f, -1, 1, 1000)
print(calc_e_trapezoid6_2)
calc_e_simpson6_2 = simpson(f, -1, 1, 1000)
print(calc_e_simpson6_2)<jupyter_output>0.0<jupyter_text>The tradional trepozoid method gives a really small number, but it's not exactly zero. The Simpson method still gives an exact answer.<jupyter_code>calc_e_simpson6_3 = simpson(f, -1, 1, 50)
print(calc_e_simpson6_3)<jupyter_output>3.552713678800501e-17<jupyter_text>But the Simpson doesn't always give us the exact answer in the case where the function is a polynomial, for example here I changed the number of subdivision to 50, and the Simpson method returns us with a really small number but it's not exactly zero. Why is it not able to give us the exact answer in this case? I don't know. --- Problem 1.3We left the most difficult part to last. We repeat the same steps as above for a new function $y = \sin{1/x}$. It's a difficult function to deal with, because the function behaves wildly when x is close to 0: we have $y = \sin{(infinity)}$ which you can think of as a function that oscillates with infinite frequency. I don't think this function even has an anti-derivative, so we don't know the exact answer we should expect. Anyway, let's give it a try. We try to plot it first.<jupyter_code>a = 0
b = 10
x = np.linspace(a, b, 10000)
y = np.sin(1/x)
plt.xlabel("x")
plt.ylabel("sin(1/x)")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()<jupyter_output>/tmp/ipykernel_373/3282203340.py:4: RuntimeWarning: divide by zero encountered in true_divide
  y = np.sin(1/x)
/tmp/ipykernel_373/3282203340.py:4: RuntimeWarning: invalid value encountered in sin
  y = np.sin(1/x)<jupyter_text>It gets really crazy around x=0. But thankfully we can still plot it, and the program even knows we are dividing by 0 at x=0. Python is really smart.<jupyter_code>def f2(x) :
    return np.sin(1/x)
# we can not put the initial value of x at zero, 
# It will produce error message: "dividing by zero!"
# So we give it a really small number
calc_e_adaptive7 = adaptive_trapezoid(f2, a=0.000001, b=10, acc=0.05, output=True)
calc_e_trapezoid7 = trapezoid(f2, 0.000001, 10, 128)
print(calc_e_trapezoid7)
calc_e_simpson7 = simpson(f2, 0.000001, 10, 128)
print(calc_e_simpson7)<jupyter_output>2.725250884384502<jupyter_text>I don't know which is more accurate, since I don't know the exact answer to this integral :P  ---Changing the range to x = [0, 2𝜋].<jupyter_code>a = 0.000001
b = np.pi*2
x = np.linspace(a, b, 10000)
y = np.sin(1/x)
plt.xlabel("x")
plt.ylabel("sin(1/x)")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()<jupyter_output><empty_output><jupyter_text>We can see more vividly the wild behavior around x=0 now as we zoomed in a little. Here we tried setting a small value for the initial value of x and avoided x = 0.<jupyter_code>calc_e_adaptive8 = adaptive_trapezoid(f2, a=0.000001, b=np.pi*2, acc=0.005, output=True)
calc_e_trapezoid8 = trapezoid(f2, 0.000001, np.pi*2, 64)
print(calc_e_trapezoid8)
calc_e_simpson8 = simpson(f2, 0.000001, np.pi*2, 64)
print(calc_e_simpson8)<jupyter_output>2.1362780653151536<jupyter_text>---Changing the range to x = [-1, 1].<jupyter_code>a = -1
b = 1
x = np.linspace(a, b, 10000)
y = np.sin(1/x)
plt.xlabel("x")
plt.ylabel("sin(1/x)")

# An integral is the area under a curve.
plt.fill_between(x, y)
plt.show()<jupyter_output><empty_output><jupyter_text>This graph actually look nice! And we can tell it is an odd function. So the integration result we expect should be 0. And surprisingly the "divide by zero" error didn't show up this time. I don't know why.<jupyter_code>calc_e_adaptive9 = adaptive_trapezoid(f2, a=-1, b=1, acc=0.005, output=True)<jupyter_output>N = 2,  Integral = 0.0
N = 2.0,  Integral = nan<jupyter_text>We can not avoid x=0 in the integration algorithm, so here is what I will do, I will divide the calculation into two parts, x = [0.000001, 1] and x = [-1, -0.000001], and add the two results together.<jupyter_code>calc_e_adaptive9_1 = adaptive_trapezoid(f2, a=0.000001, b=1, acc=0.005, output=True)
calc_e_adaptive9_2 = adaptive_trapezoid(f2, a=-1, b= -0.000001, acc=0.005, output=True)
calc_e_adaptive9 = calc_e_adaptive9_1 + calc_e_adaptive9_2
print (calc_e_adaptive9)<jupyter_output>-1.7763568394002505e-14<jupyter_text>This is a really small number close to zero. This is what we expected.<jupyter_code>calc_e_trapezoid9_1 = trapezoid(f2, 0.000001, 1, 512)
print(calc_e_trapezoid9_1)
calc_e_trapezoid9_2 = trapezoid(f2, -1, -0.000001, 512)
print(calc_e_trapezoid9_2)
calc_e_trapezoid9 = calc_e_trapezoid9_1 + calc_e_trapezoid9_2
print (calc_e_trapezoid9)<jupyter_output>-1.7541523789077473e-14<jupyter_text>The adaptive and traditional trepazoid methods now give us different values now...why? Maybe the numbers are just too small?<jupyter_code>calc_e_simpson9_1 = simpson(f2, 0.000001, 1, 512)
calc_e_simpson9_2 = simpson(f2, -1, -0.000001, 512)
calc_e_simpson9 = calc_e_simpson9_1 + calc_e_simpson9_2
print (calc_e_simpson9)<jupyter_output>-2.19824158875781e-14
{'fraction_non_alphanumeric': 0.0932025701724721, 'fraction_numerical': 0.06743320933378424, 'mean_word_length': 5.1378165213781655, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 77, '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': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '253799', 'license': 'Apache-2.0', 'n_tokens_mistral': 5962, 'n_tokens_neox': 5172, 'n_words': 1877, 'path': 'HW8_1.ipynb', 'repo_name': 'piaohanx/CP1_Calculus', 'size': '14785'}