интеграция Симпсона на питоне - PullRequest
Новая
       10

интеграция Симпсона на питоне

0 голосов
/ 14 января 2019

Я пытаюсь интегрировать численно, используя правило интеграции Симпсона для f (x) = 2x от 0 до 1, но получаю большую ошибку. Желаемое значение равно 1, но значение Python равно 1.334. Может кто-нибудь помочь мне найти решение этой проблемы? спасибо. 

import numpy as np

def f(x):
    return 2*x

def simpson(f,a,b,n):
    x = np.linspace(a,b,n)
    dx = (b-a)/n
    for i in np.arange(1,n):
        if i % 2 != 0:
            y = 4*f(x)
        elif i % 2 == 0:
            y = 2*f(x)
    return (f(a)+sum(y)+f(x)[-1])*dx/3

a = 0
b = 1
n = 1000
ans = simpson(f,a,b,n)
print(ans)

1 Ответ

0 голосов
/ 14 января 2019

Там все не так. x - это массив, каждый раз, когда вы вызываете f(x), вы оцениваете функцию по всему массиву. Поскольку n является четным и n-1 нечетным, y в последнем цикле равно 4*f(x) и из его суммы вычисляется что-то

Тогда n - количество сегментов. Количество баллов n+1. Правильная реализация 

def simpson(f,a,b,n):
    x = np.linspace(a,b,n+1)
    y = f(x)
    dx = x[1]-x[0]
    return (y[0]+4*sum(y[1::2])+2*sum(y[2:-1:2])+y[-1])*dx/3

simpson(lambda x:2*x, 0, 1, 1000)

, который затем корректно возвращает 1.000. Возможно, вы захотите добавить тест, если n четный, и увеличить его на единицу, если это не так.

Если вы действительно хотите сохранить цикл, вам нужно накапливать сумму внутри цикла. 

def simpson(f,a,b,n):
    dx = (b-a)/n;
    res = 0; 
    for i in range(1,n): res += f(a+i*dx)*(2 if i%2==0 else 4); 
    return (f(a)+f(b) + res)*dx/3;

simpson(lambda x:2*x, 0, 1, 1000)

Но циклы обычно медленнее, чем векторизованные операции, поэтому, если вы используете numpy, используйте векторизованные операции. Или просто используйте напрямую scipy.integrate.simps.

...