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

Question

Для класса числовых методов мне нужно написать программу для оценки определенного интеграла с составным правилом Симпсона. Я уже получил это далеко (см. Ниже), но мой ответ неправильный. Я тестирую программу с f (x) = x, интегрированной над 0 на 1, для которой результат должен быть равен 0,5. Я получаю 0.78746 ... и т. Д. Я знаю, что в Scipy есть правило Симпсона, но мне действительно нужно написать его сам.

Я подозреваю, что с двумя петлями что-то не так. Я пытался «для i в диапазоне (1, n, 2)» и «для i в диапазоне (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

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

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) 

будет интегрировать 2Ln (х/5) в интервале 3 и 4