Python 中的辛普森规则
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/16001157/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Simpson's rule in Python
提问by Geraldine
For a numerical methods class, I need to write a program to evaluate a definite integral with Simpson's composite rule. I already got this far (see below), but my answer is not correct. I am testing the program with f(x)=x, integrated over 0 to 1, for which the outcome should be 0.5. I get 0.78746... etc. I know there is a Simpson's rule available in Scipy, but I really need to write it myself.
对于数值方法类，我需要编写一个程序来使用辛普森复合规则计算定积分。我已经走了这么远（见下文），但我的答案不正确。我正在使用 f(x)=x 测试程序，积分超过 0 到 1，结果应该是 0.5。我得到 0.78746... 等等。我知道 Scipy 中有一个 Simpson 规则可用，但我真的需要自己编写它。
I suspect there is something wrong with the two loops. I tried "for i in range(1, n, 2)" and "for i in range(2, n-1, 2)" before, and this gave me a result of 0.41668333... etc. I also tried "x += h" and I tried "x += i*h". The first gave me 0.3954, and the second option 7.9218.
我怀疑这两个循环有问题。我之前尝试过“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)
采纳答案by unkulunkulu
You probably forget to initialize xbefore the second loop, also, starting conditions and number of iterations are off. Here is the correct way:
您可能忘记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)
Your mistakes are connected with the notion of a loop invariant. Not to get into details too much, it's generally easier to understand and debug cycles which advance at the end of a cycle, not at the beginning, here I moved the x += 2 * hline to the end, which made it easy to verify where the summation starts. In your implementation it would be necessary to assign a weird x = a - hfor the first loop only to add 2 * hto it as the first line in the loop.
你的错误与循环不变量的概念有关。不要过多地介绍细节，通常更容易理解和调试在循环结束时而不是在开始时前进的循环，在这里我将x += 2 * h行移到末尾，这样可以轻松验证求和开始的位置。在您的实现中，有必要x = a - h为第一个循环分配一个奇怪的东西，只是为了将2 * h其添加为循环中的第一行。
回答by TechnicalFox
All you need to do to make this code work is add a variable for a and b in the function bounds() and add a function in f(x) that uses the variable x. You could also implement the function and bounds directly into the simpsonsRule function if desired... Also, these are functions to be implimented into a program, not a program itself.
要使此代码工作，您需要做的就是在函数 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<<<
回答by TechnicalFox
You can use this program for calculating definite integrals by using Simpson's 1/3 rule. You can increase your accuracy by increasing the value of the variable panels.
您可以使用此程序通过使用 Simpson 的 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
For example:
例如：
def f(x,a,b):
return a * np.log(x/b)
I = integration(f,3,4,2,5)
print(I)
will integrate 2ln(x/5) within the interval 3 and 4
将在区间 3 和 4 内积分 2ln(x/5)
回答by Patryk K
There is my code (i think that is the most easy method). I done this in jupyter notebook. The easiest and most accurate code for Simpson method is 1/3.
有我的代码（我认为这是最简单的方法）。我在 jupyter notebook 中做到了这一点。辛普森方法最简单、最准确的代码是 1/3。
Explanation
解释
For standard method (a=0, h=4, b=12) and f=100-(x^2)/2
对于标准方法 (a=0, h=4, b=12) 和 f=100-(x^2)/2
We got: n= 3.0, y0 = 100.0, y1 = 92.0, y2 = 68.0, y3 = 28.0,
我们得到：n= 3.0, y0 = 100.0, y1 = 92.0, y2 = 68.0, y3 = 28.0,
So simpson method = h/3*(y0+4*y1+2*y2+y3) = 842,7 (this is not true). Using 1/3 rule we got:
所以辛普森方法 = h/3*(y0+4*y1+2*y2+y3) = 842,7（这不是真的）。使用 1/3 规则我们得到：
h = h/2= 4/2= 2 and then:
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.0, y0 = 100.0, y1 = 98.0, y2 = 92.0, y3 = 82.0, y4 = 68.0, y5 = 50.0, y6 = 28.0,
Now we calculate the integral for each step (n=3 = 3 steps):
现在我们计算每一步的积分（n=3 = 3 步）：
Integral of the first step: h/3*(y0+4*y1+y2) = 389.3333333333333
第一步的积分：h/3*(y0+4*y1+y2) = 389.3333333333333
Integral of the second step: h/3*(y2+4*y3+y4) = 325.3333333333333
第二步积分：h/3*(y2+4*y3+y4) = 325.3333333333333
Integral of the third step: h/3*(y4+4*y5+y6) = 197.33333333333331
第三步的积分：h/3*(y4+4*y5+y6) = 197.33333333333331
Sum all, and we get 912.0 AND THIS IS TRUE
总而言之，我们得到 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) )
At the beginning I added simple data entry:
一开始我添加了简单的数据输入：
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) )
回答by user1150
Example of implementing Simpson's rule for integral sinX with a = 0 and b = pi/4. And use 10 panels for the integration
对 a = 0 和 b = pi/4 的积分 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)
