数值分析之牛顿插值方法

 

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

数据点通用格式

x0x1x2xn
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−x0i=1,2,…,n
▽2yi=▽yi−▽y1xi−x1i=2,3,…,n
▽3yi=▽2yi−▽2y2xi−x2i=3,4,…,n

▽nyn=▽n−1yn−▽n−1yn−1xn−xn−1n
求常系数ai
a0=y0a1=▽y1a2=▽2y2an=▽nyn
列表分析:当n=4

x0y0    
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;">

posted @ 2018-05-19 23:59  既生喻何生亮  阅读(1569)  评论(0编辑  收藏  举报