# Simpson’s Rule (Riemann Sums with Parabolas)

int_a^b f(x) dx ~~ (h/3)((y_0+4y_1+y_2)+(y_2+4y_3+y_4)+…+(y_(n-2)+4y_(n-1)+y_n))

int_a^b f(x) dx ~~ (h/3)(y_0+4y_1+2y_2+4y_3+2y_4+…+2y_(n-2)+4y_(n-1)+y_n)

int_a^b f(x) dx ~~ sum_(j=1)^(n/2) (h/3)(f(x_(2j-2))+4f(x_(2j-1))+f(x_(2j)))

where h=(b-a)/n and the error is |E_S|<=(M(b-a)^5)/(180*n^4) and M is an upper bound for |f^((4))| on [a,b].

#### Example 1:

Evaluate int_5^7 5x dx using Simpson’s rule with n=4

x_0 through x_4 are {5, 5.5, 6, 6.5, 7} and h=2/4=1/2

Area ~~ 1/6*(5(5)+4*5(5.5)+2*5(6)+4*5(6.5)+5(7))=60

Since f'(x)=5 and f”(x)=0, then f^((4))(x)=0 and the error is zero here.

#### Example 2:

Consider int_-1^1 (x^2+5) dx with n=4

h=2/4=1/2 so x_k=-1+k/2, which is the set {-1,-1/2,0,1/2,1}

Area ~~1/6(((-1)^2+5)+4((-1/2)^2+5)+2((0)^2+5)+4((1/2)^2+5)+((1)^2+5))

Area=10 2/3

Since f^((4))(x)=0 the error here is also zero.

#### Example 3:

int_0^2 (5t^3+8t) dt with n=4

h=2/4=1/2 and x_k=0+k/2, so the x-values are {0,1/2,1,1 1/2,2}

Area ~~ 1/6((5(0)^3+8(0))+4(5(1/2)^3+8(1/2))+2(5(1)^3+8(1))+4(5(1 1/2)^3+8(1 1/2))+(5(2)^3+8(2)))

Area =36

Since f^((4))(x)=0, the error is zero here.

#### Example 4:

int_3^9 4/s^2 ds with n=4

h=6/4=3/2 and x_k=3+3k/2, so the x-values are {3, 4.5, 6, 7.5, 9}

Area ~~1/2*(4/3^2+4*4/4.5^2+2*4/6^2+4*4/7.5^2+4/9^2)

#### Example 5:

How many intervals are needed to estimate int_1^2 1/x dx to within 0.0001 square units of the actual area?

Recall the error is |E_S|<=(M(b-a)^5)/(180*n^4) and M is an upper bound for |f^((4))| on [a,b].

f(x)=x^(-1)

f'(x)=-1*x^(-2)

f”(x)=2*x^(-3)

f”(x)=-6*x^(-4)

f^((4))(x)=24*x^(-5)

f^((5))(x)=-144*x^(-6)

I found the 5th derivative to show that f^((4))(x) is decreasing, since its derivative is negative on the interval [1,2]. So, f^((4))(1)=24 is the maximum needed for the error estimation.

Plugging into |E_S|<=(M(b-a)^5)/(180*n^4):

0.0001<=(24*(1)^5)/(180*n^4)

Solving for n we get n <=((24)/(180*0.0001))^(1/4)~~6.04

A_(S_6) = 1/3*1/6*(1/1+4*1/(1 1/6)+2*1/(1 2/6)+4*1/(1 3/6)+2*1/(1 4/6)+4*1/(1 5/6)+1/2)~~0.693169

A_(S_4) = 1/3*1/4*(1/1+4*1/(1 1/4)+2*1/(1 2/4)+4*1/(1 3/4)+1/2)~~0.693254

A_(S_2) = 1/3*1/2*(1+4*1/1.5+1/2)~~.694444

ln(2)~~0.69314718

The actual error from the 6th Simpson’s sum is about 0.00002.

A_(S_20) = sum_(n=1)^10 1/(3*20)*(1/(1 + 1/20 (2 n – 2)) + 4/(1 + 1/20 (2 n – 1)) + 1/(1 + (2 n)/20)) = 5555158368718531/8014397185594800~~0.693147375

For comparison, consider 100 terms with right endpoints: sum_(k=1)^100 1/(100 (k/100 + 1))
~~0.690653430481824215252268721472608478928428119821791431288…

Or the Trapezoidal Rule with 50 intervals:
A_(T_50)=sum_(k=1)^50 1/(2*50)*(1/(1+(K-1)/50)+1/(1+k/50))~~0.693172179

#### Example 6:

Here is a picture of a parabola intersecting the reciprocal function at x=1,2,3

The parabola, y=1/6*x^2-1*x+11/6, is shown in red.
The beautiful thing about Simpson’s rule is that we don’t actually need to calculate the coefficients of any of these parabolas. The weighted sum does all of the work for us.

# A number theory story

A friend asked me to investigate how to write algorithmically generated fraction addition problems where the sum would always be reducible or irreducible. Why? It seems unfair, in a testing situation, to have one addition problem lead to a reducible sum while another is irreducible.

# Diophantine Equations and Fraction Sums

Considering the fraction addition problem a/c+b/d=(au+bv)/(LCM(c,d)), if the sum reduces then it is by a factor of the GCD(c,d). So, au+bv=n, where n is a multiple of a factor of the GCD.

# Relatively Prime Linear Combinations

If there exist x and y for which ax+by=GCD(a,b) then GCD(x,y)=1.

Let g=GCD(a,b) and x and y solutions to ax+by=g.

Since g|a and g|b, there exist integers n and m such that a=gn and b=gm.

ax+by=gnx+gmy=g

So, nx+my=1

The GCD(n,m)=1 since g is the GCD(a,b).

Let d be a divisor of x. If d also divides my, then d must divide 1, which isn’t possible. So, GCD(x,y)=1.