Правило Симпсона в Python

Вероятно, вы забыли инициализировать 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 в качестве первой строки в цикле.

Все, что вам нужно сделать, чтобы этот код заработал, - это добавить переменную для a и b в функцию bounds () и добавить функцию в 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<<<

Вы можете использовать эту программу для вычисления определенных интегралов по правилу Симпсона 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)

проинтегрирует 2ln (x / 5) в интервале 3 и 4

Вот мой код (думаю, это самый простой способ). Я сделал это в блокноте jupyter. Самый простой и точный код для метода Симпсона - 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) )

Пример реализации правила Симпсона для интеграла sinX с a = 0 и b = pi / 4. И используйте 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)