Python中的辛普森规则

Question

对于数值方法类，我需要写一个程序，用辛普森合成法则计算定积分。我已经说到这里了(见下文)，但我的答案不正确。我正在测试f(x)=x的程序，积分超过0到1，其结果应该是0.5。我得到0.78746..。我知道Scipy中有一个辛普森法则，但我真的需要自己写。

我怀疑这两个循环有问题。我以前尝试过"for i in range(1，n，2)“和"for i in range(2，n-1，2)”，结果是0.41668333...我还尝试了"x += h“和"x += i*h”。第一个给了我0.3954，第二个给了7.9218。

# Write a program to evaluate a definite integral using Simpson's rule with
# n subdivisions

from math import *
from pylab import *

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a
    for i in range(1,n/2):
        x += 2*h
        k += 4*f(x)
    for i in range(2,(n/2)-1):
        x += 2*h
        k += 2*f(x)
    return (h/3)*(f(a)+f(b)+k)

def function(x): return x

print simpson(function, 0.0, 1.0, 100)

Answer 1

您可能忘记了在第二次循环之前初始化x，而且，启动条件和迭代次数也是关闭的。下面是正确的方法：

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a + h
    for i in range(1,n/2 + 1):
        k += 4*f(x)
        x += 2*h

    x = a + 2*h
    for i in range(1,n/2):
        k += 2*f(x)
        x += 2*h
    return (h/3)*(f(a)+f(b)+k)

您的错误与循环不变量的概念有关。不需要太多的细节，通常更容易理解和调试在周期结束时推进的周期，而不是在开始时，这里我将x += 2 * h行移到了末尾，这使得验证求和从哪里开始变得很容易。在您的实现中，有必要为第一个循环分配一个奇怪的x = a - h，以便将2 * h作为循环中的第一行添加到其中。

Answer 2

要使这段代码正常工作，只需在函数bound()中为a和b添加一个变量，并在f( x )中添加一个使用变量x的函数。如果需要，还可以直接在simpsonsRule函数中实现该函数和边界。此外，这些是要实现到程序中的函数，而不是程序本身。

def simpsonsRule(n):

    """
    simpsonsRule: (int) -> float
    Parameters:
        n: integer representing the number of segments being used to 
           approximate the integral
    Pre conditions:
        Function bounds() declared that returns lower and upper bounds of integral.
        Function f(x) declared that returns the evaluated equation at point x.
        Parameters passed.
    Post conditions:
        Returns float equal to the approximate integral of f(x) from a to b
        using Simpson's rule.
    Description:
        Returns the approximation of an integral. Works as of python 3.3.2
        REQUIRES NO MODULES to be imported, especially not non standard ones.
        -Code by TechnicalFox
    """

    a,b = bounds()
    sum = float()
    sum += f(a) #evaluating first point
    sum += f(b) #evaluating last point
    width=(b-a)/(2*n) #width of segments
    oddSum = float()
    evenSum = float()
    for i in range(1,n): #evaluating all odd values of n (not first and last)
        oddSum += f(2*width*i+a)
    sum += oddSum * 2
    for i in range(1,n+1): #evaluating all even values of n (not first and last)
        evenSum += f(width*(-1+2*i)+a)
    sum += evenSum * 4
    return sum * width/3

def bounds():
    """
    Description:
        Function that returns both the upper and lower bounds of an integral.
    """
    a = #>>>INTEGER REPRESENTING LOWER BOUND OF INTEGRAL<<<
    b = #>>>INTEGER REPRESENTING UPPER BOUND OF INTEGRAL<<<
    return a,b

def f(x):
    """
    Description:
        Function that takes an x value and returns the equation being evaluated,
        with said x value.
    """
    return #>>>EQUATION USING VARIABLE X<<<

Answer 3

你可以用这个程序根据辛普森的1/3法则计算定积分。您可以通过增加可变面板的值来提高精度。

import numpy as np

def integration(integrand,lower,upper,*args):    
    panels = 100000
    limits = [lower, upper]
    h = ( limits[1] - limits[0] ) / (2 * panels)
    n = (2 * panels) + 1
    x = np.linspace(limits[0],limits[1],n)
    y = integrand(x,*args)
    #Simpson 1/3
    I = 0
    start = -2
    for looper in range(0,panels):
        start += 2
        counter = 0
        for looper in range(start, start+3):
            counter += 1
            if (counter ==1 or counter == 3):
                I += ((h/3) * y[looper])
            else:
                I += ((h/3) * 4 * y[looper])
    return I

例如：

def f(x,a,b):
    return a * np.log(x/b)
I = integration(f,3,4,2,5)
print(I)

将在区间3和4内积分2ln(x/5)