Трапецеидальное правило в Python

Question

Я пытаюсь реализовать трапецеидальное правило в Python 2.7.2. Я написал следующую функцию:

def trapezoidal(f, a, b, n): 
    h = float(b - a)/n 
    s = 0.0 
    s += h * f(a) 
    for i in range(1, n): 
     s += 2.0 * h * f(a + i*h) 
    s += h * f(b) 
    return s 

Однако, е (лямбда х: х ** 2, 5, 10, 100) возвращает 583,333 (он должен вернуться 291.667), так ясно есть что-то не так с мой сценарий. Я не могу это заметить.

Answer 1

Вы отключены в два раза. Действительно, Trapezoidal Rule, как описано в математическом классе будет использовать приращение как

s += h * (f(a + i*h) + f(a + (i-1)*h))/2.0 

(f(a + i*h) + f(a + (i-1)*h))/2.0 усредняет высоту функции в двух соседних точках сетки.

Поскольку каждая из двух соседних трапеций имеет общий край, приведенная выше формула требует оценки функции в два раза чаще, чем необходимо.

Более эффективная реализация (ближе к тому, что вы вывесили), будет сочетать общие термины из смежных итераций for-loop:

f(a + i*h)/2.0 + f(a + i*h)/2.0 = f(a + i*h) 

прибыть по адресу:

def trapezoidal(f, a, b, n): 
    h = float(b - a)/n 
    s = 0.0 
    s += f(a)/2.0 
    for i in range(1, n): 
     s += f(a + i*h) 
    s += f(b)/2.0 
    return s * h 

print(trapezoidal(lambda x:x**2, 5, 10, 100)) 

, который дает

291.66875 

Answer 2

trapezoidal rule имеет большой /2 (каждый термин равен (f(i) + f(i+1))/2, а не f(i) + f(i+1)), который вы оставили вне своего кода.

Вы использовали общую оптимизацию, которая рассматривает первую и последнюю пару специально, так что вы можете использовать 2 * f(i) вместо вычисления f(i) дважды (один раз как f(j+1) и один раз как f(i)), так что вы должны добавить / 2 к шагу петли и специальные первые и последние шаги:

s += h * f(a)/2.0 
for i in range(1, n): 
    s += 2.0 * h * f(a + i*h)/2.0 
s += h * f(b)/2.0 

Вы можете, очевидно, упростить стадию цикла, заменяя 2.0 * …/2.0 только с ….

Однако, еще проще, вы можете просто разделить все это на 2 в конце концов, ничего не меняя, но этой линии:

return s/2.0