Das deutsche Python-Forum

Code: Alles auswählen

def gauss(A):    # abkopiert von einem beispiel
    n = len(A)

    for i in range(0, n):
        # Search for maximum in this column
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i, n+1):
            tmp = A[maxRow][k]
            A[maxRow][k] = A[i][k]
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i+1, n):
            c = -A[k][i]/A[i][i]
            for j in range(i, n+1):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]

    # Solve equation Ax=b for an upper triangular matrix A
    x = [0 for i in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = A[i][n]/A[i][i]
        for k in range(i-1, -1, -1):
            A[k][n] -= A[k][i] * x[i]
    return x

def perispline(x,y):
    l = len(x)
    num_eq = l-1

    em_rows = [0]*(num_eq -1 )
    h = [x[i+1] - x[i] for i in range(l-1)]

    # first row
    r = [0.] * (num_eq -1 + 1) 
    r[0] = 2*(h[0] + h[1])
    r[1] = h[1]
    r[len(r) -2] = h[0]
    em_rows[0] = r

    # body rows
    i = 0
    while i < len(em_rows)-2:
       r = [0.] * (num_eq -1 + 1) 

       r[i+0] = h[i+1]
       r[i+1] = 2*(h[i+1]+h[i+2])
       r[i+2] = h[i+2]

       em_rows[i+1] = r
       i += 1

    # last row
    r = [0.] * (num_eq -1 + 1) 
    r[0] = h[0]
    r[len(r)-3] = h[len(h) - 2]
    r[len(r)-2] = 2 * (h[len(h) - 2] + h[len(h) - 1])
    em_rows[len(em_rows)-1] = r

    # result column
    i = 0
    while(i < len(em_rows)):
        em_rows[i][len(em_rows[i])-1] = 3. / h[i+1] * (y[i+2] - y[i+1]) - 3. / h[i+0] * (y[i+1] - y[i+0])
        i += 1

    # solve equation system for c-coeff
    gr = gauss(em_rows)

    # ich vermute das es bis hier richtig ist....

    # coeffs
    a = [0.]*num_eq
    b = [0.]*num_eq
    c = [0.]*num_eq
    d = [0.]*num_eq

    # set c and a coeffs
    i = 0
    while i < num_eq:
        if i == 0:
            c[0] = gr[-1]
        else:
            c[i] = gr[i-1]
        a[i] = y[i-1]
        i += 1

    i = 0
    while i < num_eq:        
        b[i] =  (a[(i+1)%num_eq] - a[i])/h[i]  - h[i]/3.*(c[(i+1)%num_eq] + 2*c[i])
        d[i] = (c[(i+1)%num_eq] - 2*c[i]) / (3.*h[i])
        i += 1

    def pericalc(x_value):
        #find interval
        i = 0
        while i < len(x)-1:
            if (x_value >= x[i] and x_value <= x[i+1]):
                break
            i += 1
        return a[i] + b[i]*(x_value - x[i]) + c[i]*(x_value - x[i])**2 + d[i]*(x_value - x[i])**3

    return pericalc

f = perispline([0.,1.,2.,3.,4.,5.,6.], [1.,2.,1.,2.,1.,2., 1.])
print f(0) 
print f(1) 
print f(2)
print f(3)
print f(4) 

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0,4,100)
y = []
for xx in x:
    y.append(f(xx))

plt.plot(x,y)
plt.show()

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt

def perispline(x,y):
    n = len(x)

    epsilon = [x[i+1] - x[i] for i in xrange(n-1)]
    delta = [y[i+1] - y[i] for i in xrange(n-1)]
   
    A = [[0. for i in xrange(n)] for i in xrange(n)]   # Init n x n matrix
    Y = [0.]*n                                         # Init result vector Y

    # Fill A
    for i in xrange(n-2):
        A[i][i+0] = epsilon[i]
        A[i][i+1] = 2*(epsilon[i+1] + epsilon[i])
        A[i][i+2] = epsilon[i+1]

    # Fill last two rows
    A[-2][0] = -2. * epsilon[0]
    A[-2][1] = -1. * epsilon[0]
    A[-2][-2] = 2. * epsilon[-1]
    A[-2][-1] = epsilon[-1]

    A[-1][0] = 1
    A[-1][-1] = -1

    # Fill result vector
    for i in xrange(n-2):
        Y[i] = 3. * (delta[i+1] / epsilon[i+1] - delta[i] / epsilon[i])
    Y[-2] = 3*(delta[n-2] / epsilon[n-2] - delta[0] / epsilon[0])        

    # Solve Ac=Y
    c =  np.linalg.solve(A, Y)  

    # Calculate missing coefficiants b and d
    b = [delta[i] / epsilon[i] - epsilon[i] / 3. * (2*c[i]+c[i+1])  for i in xrange(n-1)]
    d = [(c[i+1] - c[i]) / (3.*epsilon[i]) for i in xrange(n-1)]


    def spline_evaluate(x_value):
        #find interval
        i = 0
        while i < len(x)-1:
            if (x_value >= x[i] and x_value <= x[i+1]):
                break
            i += 1

        if i >= len(x)-1: 
            i = len(x)-2
        if x_value < x[0]:
            i = 0

        return y[i] + b[i]*(x_value - x[i]) + c[i]*(x_value - x[i])**2 + d[i]*(x_value - x[i])**3

    return spline_evaluate