数值分析之牛顿插值方法

 

 

数值分析之牛顿插值方法详解

数据点通用格式

x0	x1	x2	…	xn
y0	y1	y2	…	yn

算法描述

数据点有n+1个
设置常系数项n+1个，即:

[a0,a1,a2,…,an−2,an1,an]

设置系数函数Pn(x)，如下：

Pn(x)=a0+(x−x0)a1+(x−x0)(x−x2)a2+⋯+(x−x0)(x−x1)…(x−xn−1)an

利用合并同类项找出规律：
假设n=3，则P3(x)的表达式如下：

P3(x)=a0+(x−x0)a1+(x−x0)(x−x1)a2+(x−x0)(x−x1)(x−x2)a3=a0+(x−x0){a1+(x−x1)[a2+(x−x2)a3]}.

于是：假设
P0(x)=a3
P1(x)=a2+(x−x2)P0(x)
P2(x)=a1+(x−x1)P1(x)
P3(x)=a0+(x−x0)P2(x)

将其推广到一般情况：
P0(x)=an
P1(x)=an−1+(x−xn−1)P0(x)
P2(x)=an−2+(x−xn−2)P1(x)
P3(x)=an−3+(x−xn−3)P2(x)
P4(x)=an−4+(x−xn−4)P3(x)
…
Pk(x)=an−k+(x−xn−k)Pk−1(x)
…
Pn−2(x)=a2+(x−x2)Pn−3(x)
Pn−1(x)=a1+(x−x1)Pn−2(x)
Pn(x)=a0+(x−x0)Pn−1(x)
充分利用公式：yi=Pn(xi)，使得多项式穿过每一个点，利用数据，最终求得每一个常系数ai。
y0=Pn(x0)=a0
y1=Pn(x1)=a0+(x1−x0)Pn−1(x1)=a0+(x1−x0)(a1+(x1−x1)Pn−2(x1))
y2=Pn(x2)=a0+(x2−x0)a1+(x2−x0)(x2−x1)a2
…
yn−2=Pn(xn−2)=a0+(xn−2−x0)a1+⋯+(xn−2−x0)(xn−2−x1)⋯(xn−2−xn−3)an−2
yn−1=Pn(xn−1)=a0+(xn−1−x0)a1+⋯+(xn−1−x0)(xn−1−x1)⋯(xn−1−xn−2)an−1
yn=Pn(xn)=a0+(xn−x0)a1+⋯+(xn−x0)(xn−x1)⋯(xn−xn−1)an
知识点补充Introducing the divided differences
▽yi=yi−y0xi−x0， i=1,2,…,n
▽2yi=▽yi−▽y1xi−x1， i=2,3,…,n
▽3yi=▽2yi−▽2y2xi−x2， i=3,4,…,n
…
▽nyn=▽n−1yn−▽n−1yn−1xn−xn−1， n
求常系数ai，
a0=y0， a1=▽y1， a2=▽2y2， ⋯， an=▽nyn
列表分析：当n=4时

x0	y0
x1	y1	▽y1
x2	y2	▽y2	▽2y2
x3	y3	▽y3	▽2y3	▽3y3
x4	y4	▽y4	▽2y4	▽3y4	▽4y4

牛顿插值算法Python代码

import numpy as np
def coeffts(xData,yData):
    m = len(xData)
    Table = np.zeros((m,m))
    Table[:,0] = np.array(yData)
    for k in range(1,m):
        for i in range(k,m):
            Table[i,k] = (Table[i,k-1]-Table[k-1,k-1])/(xData[i]-xData[k-1])
    return np.diag(Table)


def evalPoly(a,xData,x):
    n = len(xData) - 1 # Degree of polynomial
    p = a[n]
    for k in range(1,n+1):
        p = a[n-k] + (x -xData[n-k])*p
    return p 

案例分析

import numpy as np
xData = np.array([0.15,2.3,3.15,4.85,6.25,7.95])
yData = np.array([4.79867,4.49013,4.2243,3.47313,2.66674,1.51909])
a = coeffts(xData,yData)
print(" x    yInterp")
print('-'*20)
for x in np.arange(0.0,8.1,0.5):
    y = evalPoly(a,xData,x)
    print('{:3.1f} {:9.5f}'.format(x,y))

插值结果：

<wiz_tmp_tag id="wiz-table-range-border" contenteditable="false" style="display: none;">