Ausgleich mit Polynomen

Mit Polynomen lassen sich nahezu beliebige Kurvenverläufe annähern. Alle notwendigen Daten werden im Programm selbst eingetragen.
Im selben Verzeichnis muss die Programmsammlung plot liegen.

Eine kurze Einführung in die Thematik steht als Lektion bereit.

# Ausgleich durch Polynome, piu@201119
from math import *
from sys import *
from plot import *

######## Eingaben hier ################################################

n = 3 # Grad des Polynoms
m = 6 # Anzahl Datenpunkte

# die Messwerte
x = [5.2,  4.2,  3.2,  2.3,  1.46, 0.48]     
y = [2.41, 2.87, 3.55, 4.50, 6.14, 9.35]

initKoor(0,6,1, 0,12,1) # initKoor(xMin,xMax,xScala, yMin,yMax,yScala)

########################################################################

def sq(x):
    return x*x

A = [[0 for i in range(n + 2)] for j in range(n+1)] # LGS mit rechter Seite
c = [0 for i in range(n+1)]  # die Ergebnis-Koeffizienten


# LGS aufbauen
yQ=0 # yQuer, Mittelwert
for i in range(m):
    xHochZeile = 1
    for j in range(n+1):
        xPotenz = xHochZeile 
        for k in range(n+1):
            A[j][k] = A[j][k] + xPotenz           # Summe x^(j+k) 
            xPotenz = xPotenz * x[i]
        A[j][n+1] = A[j][n+1] + xHochZeile * y[i] # Summe x^j * y
        xHochZeile = xHochZeile * x[i]
    yQ=yQ+y[i]
yQ=yQ/m # Mittelwert
print (A[1][n+1])


# LGS lösen, ohne Pivotwahl
for j in range(n+1):
    for k in range(n+1):
        f = A[k][j] / A[j][j]
        for r in range(n + 2):
            if (k != j):
                A[k][r] = A[k][r] - f * A[j][r]

# Ausgabe
print("y=")
for j in range(n+1):
    c[j] = A[j][n+1] / A[j][j]
    if (c[j] > 0):
       print("+",end="")
    print(c[j],"\t* x^",j)
print()

# Grafik
X=x[0]
dx=(x[m-1]-x[0])/1000
for i in range(1001):
    Y=0
    XP=1 # Potenz von x
    for j in range(n+1):
        Y=Y+c[j]*XP
        XP=XP*X
    Plot(X,Y,3,"red")
    X=X+dx
update()

# Bestimmtheitsmaß R²
s1=0
s2=0
for i in range(m):
    X=x[i]
    Y=y[i]
    yD=0 # yDach
    XP=1 # Potenz von x
    for j in range(n+1):
        yD=yD+c[j]*XP
        XP=XP*X
    s1=s1+sq(yD-yQ)
    s2=s2+sq(Y-yQ)
print("R2",s1/s2)

Uwe Pilz, November 2020