Python реализация составного правила Ньютона - PullRequest
Новая
       82

Python реализация составного правила Ньютона

0 голосов
/ 18 июня 2020

Я пытаюсь реализовать составное правило ньютона-кота. Проблема в том, что результат моего кода довольно сильно отличается от фактического значения.

Краткое объяснение моего кода с использованием математики c formular.

Функция weight вычисляет базис Лагранжа полином L_g для g = [0,1,2, ... n-1] для n точек данных (https://en.wikipedia.org/wiki/Lagrange_polynomial)

Функция newton_cotes вычисляет правило Котеса Ньютона путем аппроксимации функции с помощью интерполяции Лагранжа (https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas)

Функция sum_newton_cotes делит интервал интегрирования [a,b] и применяет newton_cotes к подынтервалу с k числом точки данных в этом подынтервале, а затем все суммировать. В этом коде я выбираю эквидистантные подинтервалы.

Я провожу тест на sum_newton_cotes с функцией sin(x) для [0,pi]. Для разных k и n я получаю совершенно разные результаты. Но max равен 1, и когда я увеличиваю n на 1000, результат становится очень маленьким. Любое объяснение / предложение было бы действительно полезно 

import numpy as np

def weight(x,alist):
    '''
    Input
    x: x-coordinate of point of interest
    alist: list of x-coordinates of exisiting data points

    Output
    result: Array of lagrange polynomial evaluated at x
    '''
    n = len(alist)
    result = np.zeros(n)
    for i in range(n):
        prod = 1
        for j in range(n):
            if i != j:
                prod*= (x - alist[j]) / (alist[i] - alist[j])
                result[i] = prod
    return result

def newton_cotes(a,b,n,f, dx = 0.0001):
    '''
    Input
    a,b: Integrationsintervall
    n: Anzahl von Stützstellen
    f: zu integrierende Funktion
    Output
    Approximierung des Integrals von f im [a,b]
    '''
    index = np.arange(0, n, 1)
    x_knot = [a + (u*(a-b)/n) for u in index] # n Punkte
    y_knot = [f(i) for i in x_knot]
    x_eval = np.linspace(a,b)

    # Gewicht lambda
    w = np.zeros(n)
    for i in range(n):
        w[i]=(1/(b-a))*np.sum([weight(m,x_knot)[i]*dx for m in x_eval])

    asum = 0
    for i in range(n):
        asum += y_knot[i]*w[i]

    return (b-a)*asum



def summed_newton_cotes(a,b,f,k,n):
    '''
    Input
    a, b: das Integrationsintervall
    f: die zu integrierende Funktion
    n: Anzahl der Intervalle
    k: die k-te Newton-Cotes Formel
    Output
    A: Anzahl der durchgeführten f-Auswertungen
    S: Wert der QF
    '''
    index = np.arange(0,1+((b-a)/n),1)

    # Äquidistante Intervalle
    intervall =[]
    for i in range(len(index)-1):
        intervall.append([a+index[i]*(b-a)/n,a+index[i+1]*(b-a)/n])

    result = []
    for j in range(len(intervall)):
        result.append(newton_cotes(intervall[j][0],intervall[j][1],k,f,dx=0.01))

    return np.sum(result)

def fun(x):
    return np.sin(x)

print(summed_newton_cotes(0,np.pi,fun,5,5))
...