在Python中集成

Question

我一直有意在Python中集成，但我不使用Scipy，Numpy或任何其他可以集成到python中的程序。在编码方面，我几乎是个新手，但我需要帮助整合。我已经复制了一小段代码，但我仍然需要改进它。

def LeftEndSum(startingx, endingx, numberofRectangles) : 

    width = (float(endingx) - float(startingx))/numberofRectangles 

    runningSum = 0 

    for i in range(numberofRectangles) : 
     height = f(startingx + i*width) 
     area = height * width 
     runningSum += area 
    return runningSum 

我想整合，但我想要得到的数据点的列表，然后我就可以进入图形在积分

的最后一个情节我定义的时间间隔的想法[a，b]和delta n =（在点之间的＃个盒子中的变化），其中我可以进行收敛测试以停止循环以获得点。如果I（n（旧值）+ delta（n）） - I（n（旧值））/ I（n（旧）），则收敛测试将变为

其中epsilon = 1x10^-6

，其中，如果积分值的代码打破

Answer 1

比方说，你要为y =集成X * X为0.0至10.0，我们知道答案是333.33333。下面是一些代码来做到这一点：

def y_equals_x_squared(x): 
    y = x*x 
    return y 

def LeftEndSum(startingx, endingx, numberofRectangles) : 
    width = (float(endingx) - float(startingx))/numberofRectangles 
    #print "width = " + str(width) 
    runningSum = 0 
    i = 1 
    while i <= numberofRectangles: 
     x = (endingx - startingx)/(numberofRectangles) * (i - 1) + startingx 
     height = y_equals_x_squared(x) 
     area = height * width 
     #print "i, x , height, area = " + str(i) + ", " + str(x) + ", " + str(height) + ", " + str(area) 
     runningSum += area 
     i += 1 
    return runningSum 
#----------------------------------------------------------------------------- 
startingx = 0.0 
endingx = 10.0 
# 
numberofRectangles = 3 
old_answer = LeftEndSum(startingx, endingx, numberofRectangles) 
# 
numberofRectangles = 4 
new_answer = LeftEndSum(startingx, endingx, numberofRectangles) 
# 
delta_answer = abs(new_answer - old_answer) 
# 
tolerance = 0.0001 
max_iterations = 500 
iteration_count = 0 
iterations = [] 
answers = [] 
while delta_answer > tolerance: 
    numberofRectangles += 100 
    new_answer = LeftEndSum(startingx, endingx, numberofRectangles) 
    delta_answer = abs(new_answer - old_answer) 
    old_answer = new_answer 
    iteration_count += 1 
    iterations.append(iteration_count) 
    answers.append(new_answer) 
    print "iteration_count, new_answer = " + str(iteration_count) + ", " + str(new_answer) 
    if(iteration_count > max_iterations): 
     print "reached max_iterations, breaking" 
     break 
# 
OutputFile = "Integration_Results.txt" 
with open(OutputFile, 'a') as the_file: 
    for index in range(len(answers)): 
     the_file.write(str(index) + " " + str(answers[index]) + "\n") 
# 
import matplotlib.pyplot as plt 
# 
fig, ax = plt.subplots() 
ax.plot(iterations, answers, 'r-', label = "Increasing # Rectangles") 
title_temp = "Simple Integration" 
plt.title(title_temp, fontsize=12, fontweight='bold', color='green') 
ax.legend(loc='best', ncol=1, fancybox=True, shadow=True) 
plt.xlabel('Number of Iterations') 
plt.ylabel('Answer') 
ax.grid(True) 
plt.show(block=True) 

通知我们绘制的答案与在最后迭代次数，它非常缓慢接近真实的答案，因为迭代增加数量。还有其他的集成方法比简单的矩形如梯形法则更好。当你放入一个while循环并检查一个容差时，总是要放入一个max_iterations检查，以免你陷入无限循环。

你可以在这里查看你的答案： http://www.integral-calculator.com/ 这就是我们如何知道答案是333.3333

Answer 2

它更有意义，包括f作为函数的参数。为什么选择固定配方？

def riemannSum(f,a,b,n,sample = 'L'): 
    """computes Riemann sum of f over [a,b] using n rectangles and left ('L'), right ('R') or midpoints ('M')""" 
    h = (b-a)/float(n) 
    if sample.upper() == 'L': 
     s = a #s = first sample point 
    elif sample.upper() == 'R': 
     s = a + h 
    else: 
     s = a + h/2.0 
    return h*sum(f(s+i*h) for i in range(n)) 

您可以明确定义函数，然后将它们整合：

>>> def reciprocal(x): return 1/x 

>>> riemannSum(reciprocal,1,2,100) 
0.6956534304818242 

（精确值。此外，该代码可以大大利用应用于发电机内置功能sum简化自然对数的2，这大约是0.693147）

或者，您可以使用匿名函数（lambda表达式）：

>>> riemannSum(lambda x: x**2,0,1,100,'m') 
0.333325 

或者，你已经可以在math模块中使用的功能：

>>> riemannSum(math.sin,0,math.pi,10) 
1.9835235375094546 

这些方法都不是非常准确的。更准确的是Simpson's Rule这也是很容易在Python做到：

def simpsonsRule(f,a,b,n): 
    if n%2 == 1: 
     return "Not applicable" 
    else: 
     h = (b-a)/float(n) 
     s = f(a) + sum((4 if i%2 == 1 else 2)*f(a+i*h) for i in range(1,n)) + f(b) 
     return s*h/3.0 

例如：

>>> simpsonsRule(math.sin,0,math.pi,10) 
2.0001095173150043 

这比黎曼和10米的矩形更为准确的（真正价值2）。