Блог

Для моего численного курса я должен использовать правило Симпсона для аппроксимации интеграла от трех различных функций, которые перечислены ниже…. Мой код выполняется, но он не дает мне правильного приближения, и я не знаю почему….

 import math
def F(x):
    c = x * math.log(x)
    return c
def G(x):
    g = 2 / (x**2   4)
    return g
def H(x):
    h = x**2 * math.cos(x)
    return h
def simpsonRule(f, a, b, n):
    h = (b-a)/n
    k = 0
    z = 0
    for i in range(1, (n//2) - 1):
        k  = 2*f(a (2*i*h))

    for i in range(1, (n//2)):
        z  = 4*f(a ((2*i)-1)*h)    
    return (h*(f(a)   k   z   f(b))) / 3


AA = simpsonRule(F, 1, 2, 4)
CC = simpsonRule(G, 0, 2, 6)
DD = simpsonRule(H, 0, math.pi, 6)

print(AA)
print(CC)
print(DD)
  

Я считаю, что проблема в том, что верхний предел ваших range() вызовов отключен на 1. Верхний предел никогда не достигается, поэтому он должен быть выше того, что вы хотите. Я включил дополнительный тестовый пример, sin(x) , для которого я знал правильный ответ (~ 1):

 import math

def E(x):
    return math.sin(x)

def F(x):
    return x * math.log(x)

def G(x):
    return 2 / (x**2   4)

def H(x):
    return x**2 * math.cos(x)

def simpsonRule(f, a, b, n):
    h = (b - a) / n
    k = 0
    z = 0

    for i in range(1, n // 2):
        k  = 2 * f(a   2 * i * h)

    for i in range(1, n // 2   1):
        z  = 4 * f(a   (2 * i - 1) * h)

    return h * (f(a)   k   z   f(b)) / 3.0

AA = simpsonRule(F, 1, 2, 4)
CC = simpsonRule(G, 0, 2, 6)
DD = simpsonRule(H, 0, math.pi, 6)
EE = simpsonRule(E, 0, math.pi / 2, 100)

print(AA)
print(CC)
print(DD)
print(EE)
  

ВЫВОД

 > python3 test.py
0.6363098297969493
0.7853979452340107
-6.274868388453119
1.0000000003382359
> 
  

В случае, если вы используете Python 2.x, проблема может быть вызвана типами переменных. Все переменные в правой части эксперимента h = (b-a)/n являются int , поэтому интерпретатор предполагает, что h int это также должно быть. Проблема может быть решена путем преобразования типа

 a = float(a)
b = float(b) 
  

вставлено перед этим выражением.

Упомянутое выше поведение характерно для Python 2, и оно было исправлено в Python 3, предполагая float , что в этой ситуации. Смотрите здесь больше примеров с приведением типов в Python.

Вывод без преобразования:

 0.0
0
-1.36541969235
  

Вывод с преобразованием:

 0.208501009808
0.577327327327
-1.36541969235