I have a code for Simpsons rule for definite integrals that I created myself. For some reason, this code is not outputting the correct value. My current function is y = x^3, the value of N is 4, the value of a is 0 and the value of b is 2. According to Simpsons rule, the area must then be 1/6(0^ 3(4 * (0.5 ^ 3))(2 * (1.0 ^ 3))(4 * (1.5 ^ 3))(2 3)). Calculating that by hand, the value is correctly 4. However, in my program, using the same logic, the output becomes 8 and changes based on the value of n. In my mind, the area should not change greatly just because n is changed to another positive number. It should in fact, become more precise. You can see what I mean from my code.
def trapezoid(lb, rb, r):
width = (rb-lb)/r
currentX = lb
area = currentX ** 3
yes = "4"
currentX = width
while currentX < rb:
if yes == "4":
area = (4 * (currentX ** 3))
yes = "2"
currentX = width
continue
if yes == "2":
area = (2 * (currentX ** 3))
yes = "4"
currentX = width
continue
area = ((rb) ** 3)
area *= ((rb - lb)/(3 * rb))
return area
lftBount = int(input("Enter the left bound:"))
rgtBount = int(input("Enter the right bound:"))
Repeat = int(input("How many times do you want"))
print(trapezoid(lftBount,rgtBount,Repeat))
I've run this through in my head and on paper several times, could someone point me in the direction where I went wrong?
CodePudding user response:
From https://www.intmath.com/integration/6-simpsons-rule.php it seems to me that there is only one small mistake in your code when you divide (last statement in your function) your area by the right bound:
area *= ((rb - lb)/(3 * rb))
whereas it should be divided by the number of segments, so:
area *= (rb - lb) / (3 * r)
With this modification the output produced for input values 0, 2, and 4 is indeed 4.0.
Here is a modification of your code, where I made function f(x) = x ** 3
a parameter of trapezoid
so that it can be easily changed by
another function. I also simplified a bit the body of the while
loop.
def trapezoid(lb, rb, r, f): # f = function used
width = (rb - lb) / r
currentX = lb
area = f(currentX)
multby4 = True
currentX = width
while currentX < rb:
if multby4:
area = 4 * f(currentX)
else:
area = 2 * f(currentX)
currentX = width
multby4 = not multby4 # switch bit multby4 for the next iteration
area = f(rb)
area *= (rb - lb) / (3 * r)
return area
def f(x):
return x ** 3
lftBount = int(input("Enter the left bound:"))
rgtBount = int(input("Enter the right bound:"))
Repeat = int(input("How many times do you want"))
print(trapezoid(lftBount, rgtBount, Repeat, f))