2

ニュートンの補間多項式を使用して x=8 で y を決定し、次の配列の有限分割差分を計算しようとしています。配列は

x = 0  1  2  5.5  11  13  16  18

y=  0.5  3.134  5.9  9.9  10.2  9.35  7.2  6.2

私が持っている疑似コードはhttp://imgur.com/gallery/Lm2KXxA/newにあります。答えを教えてくれる利用可能な疑似コード、アルゴリズム、またはライブラリはありますか?

また、これは matlab http://imgur.com/gallery/L9wJaEH/newでプログラムを実行する方法だと思います。Pythonでそれを行う方法がわかりません。

4

4 に答える 4

4

これがPythonコードです。関数coefは有限分割差分係数を計算し、関数Evalは特定のノードで内挿を評価します。

import numpy as np
import matplotlib.pyplot as plt

def coef(x, y):
    '''x : array of data points
       y : array of f(x)  '''
    x.astype(float)
    y.astype(float)
    n = len(x)
    a = []
    for i in range(n):
        a.append(y[i])

    for j in range(1, n):

        for i in range(n-1, j-1, -1):
            a[i] = float(a[i]-a[i-1])/float(x[i]-x[i-j])

    return np.array(a) # return an array of coefficient

def Eval(a, x, r):

     ''' a : array returned by function coef()
        x : array of data points
        r : the node to interpolate at  '''

    x.astype(float)
    n = len( a ) - 1
    temp = a[n] + (r - x[n])
    for i in range( n - 1, -1, -1 ):
        temp = temp * ( r - x[i] ) + a[i]
    return temp # return the y_value interpolation
于 2014-11-28T14:04:08.093 に答える
0

@Ledruid によって提案されたソリューションが最適です。2 次元配列は必要ありません。分割された差分ツリーは、2D 配列のようなものです。より原始的なアプローチ (そこから @Leudruid のアルゴリズムを取得できます) は、ツリーのノードに対応する「行列 $F_{ij}$ を考えることです。このようにして、アルゴリズムは次のようになります。

import numpy as np
from numpy import *

def coeficiente(x,y) :
    ''' x: absisas x_i 
        y : ordenadas f(x_i)'''

    x.astype(float)
    y.astype(float)
    n = len(x)
    F = zeros((n,n), dtype=float)
    b = zeros(n) 
    for i in range(0,n):
        F[i,0]=y[i]



    for j in range(1, n):
        for i in range(j,n):
            F[i,j] = float(F[i,j-1]-F[i-1,j-1])/float(x[i]-x[i-j])


    for i in range(0,n):
        b[i] = F[i,i]

    return np.array(b) # return an array of coefficient

def Eval(a, x, r):

    '''  a : retorno de la funcion coeficiente() 
         x : abcisas x_i
         r : abcisa a interpolar
    '''

    x.astype(float)
    n = len( a ) - 1
    temp = a[n]
    for i in range( n - 1, -1, -1 ):
        temp = temp * ( r - x[i] ) + a[i]
    return temp # return the y_value interpolation
于 2017-09-10T21:40:51.170 に答える