牛顿插值法的Python程序

牛顿插值

差商及其性质

牛顿插值公式

#牛顿插值，差商
from sympy import *

x = Symbol('x')
x_list = [1, 3, 2]
y_list = [1, 2, -1]

n = len(x_list)
df_list = [y_list[0]]
for k in range(1,n):
    df = 0
    for i in range(k+1):
        X = 1
        for j in range(k+1):
            if i == j:
                pass
            else:
                X = X * (x_list[i]-x_list[j])
        df += y_list[i] / X
    df_list.append(df)
print(df_list)
P = df_list[0]
n = len(df_list)
for i in range(n-1):
    X = 1
    for j in range(i+1):
        X = X * (x - x_list[j])
    P += df_list[i+1] * X
print(expand(P))
print(P.subs(x,1.5))

结果：

[1, 0.5, 2.5]
2.5*x**2 - 9.5*x + 8.0
-0.625000000000000

#牛顿插值
from sympy import *

def f(a):
    f = sinh(a)
    return f
x_list = [0, 0.20, 0.30, 0.50]
y_list = []
for x in x_list:
    y = f(x)
    y_list.append(y)

n = len(x_list)
df_list = [y_list[0]]
for k in range(1,n):
    df = 0
    for i in range(k+1):
        X = 1
        for j in range(k+1):
            if i == j:
                pass
            else:
                X = X * (x_list[i]-x_list[j])
        df += y_list[i] / X
    df_list.append(df)
print(df_list)
P = df_list[0]
n = len(df_list)
x = Symbol('x')
a = 0.23
for i in range(n-1):
    X = 1
    for j in range(i+1):
        X = X * (x - x_list[j])
    P += df_list[i+1] * X
    print(f'{i+1}次牛顿差值多项式：',P,';',f'f({a})={P.subs(x,a)}')

结果：

[0, 1.00668001270547, 0.0838763211833875, 0.172461698783472]
1次牛顿差值多项式： 1.00668001270547*x ; f(0.23)=0.231536402922258
2次牛顿差值多项式： 0.0838763211833875*x*(x - 0.2) + 1.00668001270547*x ; f(0.23)=0.232115149538423
3次牛顿差值多项式： 0.172461698783472*x*(x - 0.3)*(x - 0.2) + 0.0838763211833875*x*(x - 0.2) + 1.00668001270547*x ; f(0.23)=0.232031850537911

差商和导数

from sympy import *

x_list = [0, 1, 2, 3, 4, 5]
y_list = [-7, -4, 5, 26, 65, 128]
DF = []
DF.append(y_list)
i = 0
while True:
    n = len(DF[i])
    i += 1
    df_list = []
    num = 0
    for k in range(n-1):
        df = (DF[i-1][k+1]-DF[i-1][k]) / i
        if df == 0:
            num += 1
        df_list.append(df)
    if len(df_list) == 1:
        break
    if len(df_list) == num:
        break
    DF.append(df_list)
print(DF)

F = DF[0][0]
x = Symbol('x')
for i in range(1,len(DF)):
    X = 1
    for j in range(i):
        X = X * (x - x_list[j])
    F += DF[i][0] * X
print(f'{F} = {expand(F)}')

结果：

[[-7, -4, 5, 26, 65, 128], [3.0, 9.0, 21.0, 39.0, 63.0], [3.0, 6.0, 9.0, 12.0], [1.0, 1.0, 1.0]]
1.0*x*(x - 2)*(x - 1) + 3.0*x*(x - 1) + 3.0*x - 7 = 1.0*x**3 + 2.0*x - 7

差分

等距节点牛顿插值公式

#差分，向前差分
from sympy import *
x = 0.57891
x_list = [0.4, 0.5, 0.6, 0.7]
def f(a):
    f = sin(a)
    return f
y_list = []
for x_a in x_list:
    y = f(x_a)
    y_list.append(y)
DF = [y_list]
n = len(x_list)
for i in range(n-1):
    df_list = []
    for j in range(len(DF[-1])-1):
        df = DF[i][j+1] - DF[i][j]
        df_list.append(df)
    DF.append(df_list)
print(DF)
x_list.append(x)
x_list.sort()
n_x = x_list.index(x)
h =0.1
t = (x - x_list[n_x-1]) / h
P = DF[0][n_x-1]
for i in range(1,len(DF)-1):
    X, Y= 1, 1
    for j in range(i):
        X = X * (t - j)
        Y = Y * (j+1)
    P += (DF[i][n_x-1]) * X / Y
print(P)

结果：

[[0.389418342308651, 0.479425538604203, 0.564642473395035, 0.644217687237691], [0.0900071962955525, 0.0852169347908324, 0.0795752138426556], [-0.00479026150472012, -0.00564172094817672], [-0.000851459443456604]]
0.547139672804571

#差分，向后差分
from sympy import *
x = 0.57891
x_list = [0.4, 0.5, 0.6, 0.7]
def f(a):
    f = sin(a)
    return f
y_list = []
for x_a in x_list:
    y = f(x_a)
    y_list.append(y)
DF = [y_list]
n = len(x_list)
for i in range(n-1):
    df_list = []
    for j in range(len(DF[-1])-1):
        df = DF[i][j+1] - DF[i][j]
        df_list.append(df)
    DF.append(df_list)
print(DF)
x_list.append(x)
x_list.sort()
n_x = x_list.index(x)
h =0.1
t = (x - x_list[n_x+1]) / h
P = DF[0][n_x]
for i in range(1,len(DF)-1):
    X, Y= 1, 1
    for j in range(i):
        X = X * (t - j)
        Y = Y * (j+1)
    P += (DF[i][n_x-i]) * X / Y
print(P)

结果：

[[0.389418342308651, 0.479425538604203, 0.564642473395035, 0.644217687237691], [0.0900071962955525, 0.0852169347908324, 0.0795752138426556], [-0.00479026150472012, -0.00564172094817672], [-0.000851459443456604]]
0.546058556206317