7

对于数值方法类，我需要编写一个程序来使用辛普森的复合规则评估定积分。我已经做到了这一点（见下文），但我的答案不正确。我正在用 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)
4

6 回答 6

9

您可能忘记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其添加为循环中的第一行。

于 2013-04-14T16:27:06.060 回答
2

要使此代码正常工作，您需要做的就是在函数 bounds() 中为 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<<<
于 2013-10-16T15:50:50.987 回答
1

您可以使用此程序通过使用辛普森的 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)

于 2016-09-13T15:54:50.103 回答
1

有我的代码（我认为这是最简单的方法）。我在 jupyter notebook 中完成了这个。Simpson 方法最简单、最准确的代码是 1/3。

解释

对于标准方法 (a=0, h=4, b=12) 和 f=100-(x^2)/2

我们得到：n= 3.0, y0 = 100.0, y1 = 92.0, y2 = 68.0, y3 = 28.0,

所以辛普森方法 = h/3*(y0+4*y1+2*y2+y3) = 842,7（这不是真的）。使用 1/3 规则，我们得到：

h = h/2= 4/2= 2 然后：

n= 3.0, y0 = 100.0, y1 = 98.0, y2 = 92.0, y3 = 82.0, y4 = 68.0, y5 = 50.0, y6 = 28.0,

现在我们计算每一步的积分（n=3 = 3 步）：

第一步的积分：h/3*(y0+4*y1+y2) = 389.3333333333333

第二步的积分：h/3*(y2+4*y3+y4) = 325.3333333333333

第三步的积分：h/3*(y4+4*y5+y6) = 197.33333333333331

总和，我们得到 912.0，这是真的

x=0
b=12
h=4
x=float(x)
h=float(h)
b=float(b)
a=float(x)
def fun(x): 
    return 100-(x**2)/2
h=h/2
l=0  #just numeration
print('n=',(b-x)/(h*2))
n=int((b-a)/h+1)
y = []   #tablica/lista wszystkich y / list of all "y"
yf = []
for i in range(n):
    f=fun(x)
    print('y%s' %(l),'=',f)
    y.append(f)
    l+=1
    x+=h
print(y,'\n')
n1=int(((b-a)/h)/2)  
l=1
for i in range(n1):
    nf=(h/3*(y[+0]+4*y[+1]+y[+2]))
    y=y[2:]  # with every step, we deleting 2 first "y" and we move 2 spaces to the right, i.e. y2+4*y3+y4
    print('Całka dla kroku/Integral for a step',l,'=',nf)
    yf.append(nf)
    l+=1
print('\nWynik całki/Result of the integral =', sum(yf) )

一开始我添加了简单的数据输入：

d=None
while(True):
    print("Enter your own data or enter the word "test" for ready data.\n")
    x=input ('Enter the beginning of the interval (a): ') 
    if x == 'test':
        x=0
        h=4  #krok (Δx)
        b=12 #granica np. 0>12  
        #w=(20*x)-(x**2)  lub   (1+x**3)**(1/2)
        break
    h=input ('Enter the size of the integration step (h): ')
    if h == 'test':
        x=0
        h=4 
        b=12 
        break
    b=input ('Enter the end of the range (b): ')
    if b == 'test':
        x=0
        h=4  
        b=12 
        break 
    d=input ('Give the function pattern: ')
    if d == 'test':
        x=0
        h=4  
        b=12
        break
    elif d != -9999.9:
        break

x=float(x)
h=float(h)
b=float(b)
a=float(x)

if d == None or d == 'test':
    def fun(x): 
        return 100-(x**2)/2 #(20*x)-(x**2)
else:
    def fun(x): 
        w = eval(d)
        return  w
h=h/2
l=0  #just numeration
print('n=',(b-x)/(h*2))
n=int((b-a)/h+1)
y = []   #tablica/lista wszystkich y / list of all "y"
yf = []
for i in range(n):
    f=fun(x)
    print('y%s' %(l),'=',f)
    y.append(f)
    l+=1
    x+=h
print(y,'\n')
n1=int(((b-a)/h)/2)  
l=1
for i in range(n1):
    nf=(h/3*(y[+0]+4*y[+1]+y[+2]))
    y=y[2:]
    print('Całka dla kroku/Integral for a step',l,'=',nf)
    yf.append(nf)
    l+=1
print('\nWynik całki/Result of the integral =', sum(yf) )
于 2018-10-05T15:43:32.153 回答
1 
def simps(f, a, b, N):   # N must be an odd integer
    """ define simpson method, a and b are the lower and upper limits of
    the interval, N is number of points, dx is the slice
        """
    integ = 0
    dx = float((b - a) / N)

    for i in range(1,N-1,2):
        integ += f((a+(i-1)*dx)) + 4*f((a+i*dx)) + f((a+(i+1)*dx))
        integral = dx/3.0 * integ
    
    # if number of points is even, then error araise
    if (N % 2) == 0:
      raise ValueError("N must be an odd integer.")
        

    return integral


def f(x):
    return x**2


integrate = simps(f,0,1,99)
print(integrate)
于 2020-07-18T21:11:49.593 回答
0

用 a = 0 和 b = pi/4 实现 Simpson 积分 sinX 规则的示例。并使用 10 个面板进行集成

def simpson(f, a, b, n):
  x = np.linspace(a, b, n+1)
  w = 2*np.ones(n+1); w[0] = 1.0; w[-1] = 1;
  for i in range(len(w)):
      if i % 2 == 1:
          w[i] = 4
  width = x[1] - x[0]
  area = 0.333 * width * np.sum( w*f(x))
  return area

f = lambda x: np.sin(x)
a = 0.0; b = np.pi/4

areaSim = simpson(f, a, b, 10)
print(areaSim)
于 2019-04-15T12:52:30.583 回答