Erstellen eines Interpolationspolynoms nach Newton

[Hallo810]

Hallo,
ich bin gerade am Programmieren eines Interpolationspolynoms nach Newton.
Mein Code funktioniert für die auskommentierten Werte in der Main-Methode. Jedoch leider nicht für die Werte, die ich versuche über eine Excel-Liste reinzugeben. Ich erhalte zwar den Graphen der Werte aus der Excel-Liste, aber nicht das Interpolationspolynom an den Stützstellen. Außerdem bekomme ich folgende Fehlermeldung:

e:/Dorothea/Python Skripte/test2.py:34: RuntimeWarning: overflow encountered in double_scalars
counter = counter * (X - X[k-1])
e:/Dorothea/Python Skripte/test2.py:75: RuntimeWarning: invalid value encountered in double_scalars
C = (Y - sum) / A
e:/Dorothea/Python Skripte/test2.py:74: RuntimeWarning: invalid value encountered in double_scalars
sum = sum + A[j] * C[j]

Mein Code lautet folgendermaßen:
import pandas as pd
import numpy as np
import numpy.linalg
import matplotlib.pyplot as plt

#Datensatz importieren
dataset = pd.read_excel('E:\Maus1_gesund.xlsx')
values = dataset.values

valuesX = np.zeros(len(values))
for i in range(len(values)):
valuesX = values[0]

valuesY = np.zeros(len(values))
for j in range(len(values)):
valuesY[j] = values[j][1]

#print(valuesY[:])

def newton_matrix(X):
# """Setup the matrix of the LSE which is used to determine the coefficients
# of the Newton-basis. X are the x-coordinates of the nodes which are
# used for interpolation."""
# Matrix zweidimensional mit Nullen fuellen und 1. Spalte mit 1
result = np.zeros((len(X), len(X)))
for a in range(0, len(result)):
result[a][0] = 1
# zeilenmaessiges Befuellen der Matrix mit den entsprechenden Polynomen
for i in range(1, len(result)):
for j in range(1, i+1):
counter = 1
for k in range(1, j+1):
counter = counter * (X - X[k-1])
result[j] = counter

return result



def newton_polynomial(C, X):
# """Take coefficients and interpolation point x-coordinates of the
# Newton-polynomial and determine the corresponding interpolation polynomial."""
# Bedingung muss erfuellt sein, ansonsten Abbruch der Methode
assert len(C) == len(X)

# Anlegen des Newton-Polynom; Methode poly1d erstellt ein Polynom mithilfe von True gibt er die faktorisierte
# Form an
result = np.poly1d([])

# Erstellen des Interpolationspolynoms nach Newton; np.poly1d gibt fuer eine leere Liste 1 heraus
for i in range(0, len(C)):
X_new = []
for j in range(0, i):
X_new.append(X[j])
X_poly = np.poly1d(X_new, True)
result = result + C[i] * X_poly

return result


def interpolating_polynomial(X,Y):
#"""Determine the interpolating polynomial for the given NumPy arrays of x and y coordinates."""
# Pruefen, ob Bedingung erfuellt wird, ansonsten Abbruch der Methode
assert len(X) == len(Y)
# Erstellen der Matrix A
A = newton_matrix(X)
# Erstellen des Vektors C
C = np.zeros(len(A))
# Koeffizienten des Newton-Polynoms bestimmen und als Polynom darstellen
for i in range(len(A)):
sum = 0
for j in range(len(A)):
sum = sum + A[i][j] * C[j]
C[i] = (Y[i] - sum) / A[i][i]
result = newton_polynomial(C, X)
return result


def interpolation_plot(X,Y):
p = interpolating_polynomial(X, Y)
px = np.arange(min(X)-0.1, max(X)+0.11, 0.01)
plt.grid(True)
plt.plot(X, Y, "o")
plt.plot(px, p(px))
plt.show()


def main():
# newton_matrix(np.array([-2., 0., 2.]))

# X = np.array([0, 1, 2,3])
# Y = np.array([-2.,3.,1.,2.])

X = valuesX
Y = valuesY
interpolation_plot(X, Y)

if __name__ == "__main__": main()

Meine Werte in der Excel-Liste besteht aus 215 x und y-Werten

Hat jemand eine Idee?


Vielen Dank schonmal.

[__blackjack__]

@Hallo810: Dann müsstest Du den Warnungen mal nachgehen warum denn da ein Überlauf passiert und was die ungültigen Werte sind.

Anmerkungen zum Quelltext: Ganz allgemein kann man sagen das Numpy nicht genutzt wird und der Code der Numpy versucht zu umgehen ziemlich unpythonisch aussieht.

`numpy.linalg` wird importiert aber nicht verwendet.

Auf Modulebene sollte nur Code stehen der Konstanten, Funktionen, und Klassen definiert. Das Hauptprogramm steht üblicherweise in einer Funktion die `main()` heisst. Die Funktion hast Du ja sogar und trotzdem definierst Du globale Variablen.

Namen werden in Python klein_mit_unterstrichen geschrieben. Ausnahmen sind Konstanten (KOMPLETT_GROSS) und Klassen (MixedCase).

Kommentare sollen dem Leser einen Mehrwert über den Code geben. Faustregel: Kommentare beschreiben nicht *was* der Code macht, denn das steht da bereits als Code, sondern warum er das macht. Sofern das nicht offensichtlich ist. Offensichtlich ist in aller Regel auch was in der Dokumentation von Python und den verwendeten Bibliotheken steht. So etwas wie ``# Erstellen der Matrix A`` ist mehr als überflüssig wenn da bloss eine Zeile folgt in der `A` das Ergebnis einer Funktion `newton_matrix()` zugewiesen wird. Der Kommentar enthält keine Information die nicht schon in dieser einen Zeile Code stecken würde.

``for i in range(len(sequence))`` um dann mit `i` auf die Elemente zuzugreifen ist in Python ein „anti pattern“. Man kann in Python direkt über die Elemente von Sequenzwerten iterieren, ohne den unnötigen Umweg über einen Laufindex. Auch das anlegen von ”leeren” Datenstrukturen bei denen man danach mit einer Schleife die Dummywerte durch die tatsächlichen Werte macht man in Python nicht.

Zudem sind Python-Schleifen bei Numpy ein Warnzeichen das man etwas falsch macht, denn Numpy verwendet man ja gerade weil man damit die expliziten Schleifen loswerden will und durch die vektorisierten Operationen auf Numpy-Arrays ersetzt.

Das fängt schon beim Laden der Daten und Element für Element kopieren in ein mit Nullen gefülltes Array an. An der Stelle würde man einfach die komplette Spalte aus dem Ausgangsarray ”slicen”. Oder man transponiert das Array und löst das dann durch eine simple Zuweisung an die beiden Variablen.

`min()` und `max()` zu verwenden um bei einem *sortierten* Array das erste und das letzte Element zu bestimmen ist ineffizient. Man kann da doch *direkt* drauf zugreifen.

Warum sind die Docstrings alle auskommentiert?

`sum()` ist der Name einer eingebauten Funktion, den sollte man nicht an andere Werte binden. Insbesondere an einer Stelle wo man die `sum()`-Funktion sogar gebrauchen könnte.

Wenn eine Vorbedingung erfüllt sein *muss*, dann ist `assert` das falsche Mittel. Das ist dazu da um Bedingungen zu prüfen die garantiert wahr sind. Und es muss nicht ausgeführt werden. Wenn man also eine Funktion abbrechen will, dann muss man eine Ausnahme auslösen.

In `newton_matrix()` läuft `k` von 1 bis j+1 von `k` wird beim verwenden dann aber 1 abgezogen. Warum also nicht von 0 bis j?

Zwischenstand bei dem einiges unpythonisches entfernt ist, aber Numpy sicher noch besser genutzt werden kann:

Code: Alles auswählen

#!/usr/bin/env python3
from functools import reduce
from operator import mul

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


def product(values):
    return reduce(mul, values, 1)


def newton_matrix(X):
    """
    Setup the matrix of the LSE which is used to determine the coefficients of
    the Newton-basis.  X are the x-coordinates of the nodes which are used for
    interpolation.
    """
    result = np.zeros((len(X), len(X)))
    result[:, 0] = 1
    for i, (result_row, x_row) in enumerate(zip(result[1:], X[1:]), 1):
        for j in range(1, i + 1):
            result_row[j] = product(x_row - X[k] for k in range(j))

    return result


def newton_polynomial(C, X):
    """
    Take coefficients and interpolation point x-coordinates of the
    Newton-polynomial and determine the corresponding interpolation polynomial.
    """
    assert len(C) == len(X)
    return sum(
        (c * np.poly1d(X[:i], True) for i, c in enumerate(C)), np.poly1d([])
    )


def interpolating_polynomial(X, Y):
    """
    Determine the interpolating polynomial for the given Numpy arrays of x and y
    coordinates.
    """
    assert len(X) == len(Y)
    A = newton_matrix(X)
    #
    # Koeffizienten des Newton-Polynoms bestimmen und als Polynom darstellen.
    #
    C = np.zeros(len(A))
    for i, (a_row, y_row) in enumerate(zip(A, Y)):
        C[i] = (y_row - sum(a_row * C)) / a_row[i]

    return newton_polynomial(C, X)


def interpolation_plot(X, Y):
    p = interpolating_polynomial(X, Y)
    xs = np.arange(X[0] - 0.1, X[-1] + 0.11, 0.01)
    plt.grid(True)
    plt.plot(X, Y, "o")
    plt.plot(xs, p(xs))
    plt.show()


def main():
    # X = np.array([0, 1, 2, 3])
    # Y = np.array([-2.0, 3.0, 1.0, 2.0])

    X, Y = pd.read_excel(R"E:\Maus1_gesund.xlsx").values.transpose()
    interpolation_plot(X, Y)


if __name__ == "__main__":
    main()