I am trying to integrate a polynomial function via the Trapezoidal Method (I can change to a more accurate method later). My code isn't perfect, and I'd like to understand exactly why it doesn't work. One problem I have is that the while loop does not end. My code thus far is as follows.
def Integrate_Trapezoidal(x_LoBound,x_HiBound,N):
"""
INPUT :
x_LoBound -- lower bound of integral
x_HiBound -- upper bound of integral
N -- number of slices (N --> inf ==> integral)
OUTPUT :
-- approximate value of integral
"""
## CREATE ALPHABET
alphabet = [chr(i) for i in range(ord('a'),ord('z')+1)]
## alphabet = ['a','b','c',...,'z'] ##
## WOULD LOVE TO TRY FLOATING INPUTS VIA ARRAY COMPREHENSION
a = float(input("What is the coefficient of the lowest order term: "))
CoeffList = []
CoeffNumList = []
LengthCoeffList = [] ## [1,2,3,...,max] where max = coefficient of highest-order term
for letter in alphabet:
AddOne = int(1)
AddOne += int(1)
for i in range(int(1),int(AddOne)):
letter = alphabet[int(i)]
while letter in alphabet:
CoeffList.append(letter)
LengthCoeffList.append(len(CoeffList))
# alphabet[i]
# i = i + 1
letter = float(input("What is the coefficient of the next-order term: ")) ## GO FROM a = ___ TO b = ___ TO c = ___ ...
CoeffNumList.append(letter)
if float(input("What is the coefficient of the next-order term: ")) == '0':
print("Type 'Y for YES and 'N' for NO")
YESorNO = str(input("Is that the last term of the polynomial: "))
endterm = YESorNO[-1] ## look at last character of string
if endterm == 'N' or endterm == 'n' or endterm == 'no' or endterm == 'NO' or endterm == 'No':
pass
elif endterm == 'Y' or endterm == 'y' or endterm == 'YES' or endterm == 'yes' or endterm == 'Yes':
break
def f(x):
"""
INPUT :
x -- variable of function
EX: x = x_LoBound OR x = x_HiBound
OUTPUT :
function -- f(x) = a x^0 + b x^1 + ...
EX: f(x_LoBound) OR f(x_HiBound)
"""
for expval in LengthCoeffList and CoeffNum in CoeffNumList:
# function = 0
function += CoeffNum * x**expval
return function
letter = alphabet[int(i+1)] ## GO FROM a TO b TO c ...
## TRAPEZOIDAL RULE
# def f(x):
# return x**4 - 2*x + 1
ht = (x_HiBound - x_LoBound) / N
ss = 0.5 * f(x_LoBound) + 0.5 * f(x_HiBound)
for num in range(1,N):
ss += f(x_LoBound + num*ht)
return ht*ss
checkanswer = Integrate_Trapezoidal(0,2,10)
print(checkanswer)
I've had a go at looking over your code and found something that I think works, checking against a couple of college handouts I downloaded. As you have said in the comments, there were a lot of extra lists which aren't necessary, so I've cut back the code a lot there.
In particular, if the presumption if that each coefficient is added in sequence from lowest to highest order, and 0 is added for any that aren't there, all you need is the number of the element in the list to know the power of x.
I also moved the definition of f() to create the helper function solve_point() which works the same I think. In particular, sum and enumerate are built in, with enumerate iterating through coeff_list and also returning a count to give the power (0 upwards).
get_coefficients() was from your old Integrate_Trapezoidal() but more focused on just one thing - which is why it then returns CoeffList to be finally processed at the end.
def solve_point(x, coeff_list):
return sum(coeff * x**e for e, coeff in enumerate(coeff_list))
def get_coefficients():
CoeffList = []
while True:
# GO FROM a = ___ TO b = ___ TO c = ___ ...
coeff = float(input("What is the coefficient of the next-order term: "))
CoeffList.append(coeff)
if coeff == 0:
YESorNO = raw_input("Is that the last term of the polynomial: [Y/N] ")
if YESorNO.upper() == 'Y':
return CoeffList[:-1]
lo, hi, n = 0, 2, 6
coeff_list = get_coefficients()
ht = (hi - lo) / float(n)
ss = 0.5 * solve_point(lo, coeff_list) + 0.5 * solve_point(hi, coeff_list)
for num in range(1,n):
ss += solve_point(lo + num*ht, coeff_list)
checkanswer = ht*ss
print(checkanswer)
I think its right - I have done a couple of checks. Hopefully it may be of help for your rewrite! If you have any examples that don't work, it would be good to know, or any errors you can see...