《计算方法》---拉格朗日插值法和牛顿插值法（python实现）

本文是相关插值法的代码实现，通过以下的例子

**给出下列数据表，用拉格朗日插值法和牛顿插值法编制程序，计算x取0.5,0.7和0.85时的函数近似值
i 0 1 2 3 4
xi 0.40 0.55 0.65 0.80 0.90
y 0.41075 0.57815 0.69675 0.88811 1.02652 **
方法一：拉格朗日插值法

x = [0.85]        #需要计算函数近似值的x
x0 = [0.40,0.55,0.65,0.80,0.90]
y = [0.41075,0.57815,0.69675,0.88811,1.02652]
m = len(x0)      #给的函数值的个数
print(x)
print(x0)
print(y)
print(n)

output:
[0.85]
[0.4, 0.55, 0.65, 0.8, 0.9]
[0.41075, 0.57815, 0.69675, 0.88811, 1.02652]
4

l = 1
L = 0
for i in range(0,m):
    for n in range(0,m):
        if n != i:
            temp = (x[0] - x0[n])/(x0[i] - x0[n])
            l = l * temp
    L = y[i]*l + L
    l = 1            #就是这一小步，导致结果偏差！！每次循环记得l重新赋值
print(L)

方法二：牛顿插值法

import numpy as np
a = np.array([[0.41075,0.57815,0.69675,0.88811,1.02652],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0]])
x = np.array([0.40,0.55,0.65,0.80,0.90])
x1 = 0.5                    #需要计算函数近似值的x
print(a)
print(x)

output:
[[ 0.41075 0.57815 0.69675 0.88811 1.02652]
[ 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. ]
[ 0. 0. 0. 0. 0. ]]
[ 0.4 0.55 0.65 0.8 0.9 ]

#计算差商表（通过多维数组***）
for i in range(1,5):       #从第二行循环到最后一行
    Y = 5 - i              #重要
    for j in range(0,Y):  #从第1列循环到Y-1列
        a[i][j] = (a[i-1][j+1] - a[i-1][j])/(x[j+1+i-1] - x[j])
print(a)

output:
[[ 0.41075 0.57815 0.69675 0.88811 1.02652 ]
[ 1.116 1.186 1.27573333 1.3841 0. ]
[ 0.28 0.35893333 0.43346667 0. 0. ]
[ 0.19733333 0.21295238 0. 0. 0. ]
[ 0.0312381 0. 0. 0. 0. ]]

#根据牛顿插值法公式计算函数近似值
temp = 1                    #通过临时变量temp来进行累加
y = 0
for i in range(0,4):
    temp = temp * (x1 - x[i])
    y1 = temp * a[i+1][0]
    y = y + y1
print(y+a[0][0])

问题汇集：
1.

#第一次出现的问题
temp = 1                    
y = 0
for i in range(0,4):
    y = a[0][0]                    #问题所在！这样，每一次进行循环y都会重新复制为a[0][0],不符合我意思
    temp = temp * (x1 - x[i])
    y1 = temp * a[i+1][0]
    y = y + y1                  #问题所在
print(y)

unsupported operand type(s) for -: ‘list’ and ‘float’

x = [0.40,0.55,0.65,0.80,0.90]  #改为：x = np.array([0.40,0.55,0.65,0.80,0.90])
temp = temp * (x - x[i])    

‘int’ object is not subscriptable

for i in range(1,5):
    x = 5 - i             #报错原因：不能用x，之前x用过做数组了
    for j in range(0,x):
        a[i][j] = (a[i-1][j+1] - a[i-1][j])/(x[1] - x[0])
print(a)