东南大学研究生-数值分析上机题（2023）Python 4 多项式插值与函数最佳逼近

3次样条插值函数

4.1 题目

(1) 编写求第一型3次样条插值函数的通用程序；
(2) 已知汽车门曲线型值点的数据如下：

i	0	1	2	3	4	5	6	7	8	9	10
xi	0	1	2	3	4	5	6	7	8	9	10
yi	2.51	3.30	4.04	4.70	5.22	5.54	5.78	5.40	5.57	5.7	5.80

端点条件为$y_0^{\prime }=0.8$，$y_{10}^{\prime }=0.2$，用所编程序求车门的3次样条插值函数$S(x)$，并打印出$S(i+0.5)$，$i=0,1,\cdots ,9$.

4.2 Python源程序

4.2.1 第一型3次样条插值函数

from  sympy  import *  
  
n = input("请输入插值点个数：")  
n = int(n)  
  
x_n = input("请输入插值节点的横坐标，以空格分隔：").split()  
y_n = input("请输入插值节点的函数值，以空格分隔：").split()  
  
x_n = [float(x) for x in x_n]  
y_n = [float(y) for y in y_n]  
  
dy_0 = input("请输入左边界条件:y'(x0)=")  
dy_n = input("请输入右边界条件:y'(xn)=")  
  
dy_0 = float(dy_0)  
dy_n = float(dy_n)  
  
f1, f2 = [0 for x in range(n-1)], [0 for x in range(n-2)]  
# 求d  
for i in range(n-1):  
    f1[i] = (y_n[i+1] - y_n[i]) / (x_n[i+1] - x_n[i])  
for i in range(n-2):  
    f2[i] = (f1[i+1] - f1[i]) / (x_n[i+2] - x_n[i])  
f2.insert(0, (f1[0] - dy_0)/(x_n[1] - x_n[0]))  
f2.append((dy_n - f1[n-2])/(x_n[n-1] - x_n[n-2]))  
d = [x * 6 for x in f2]  
# 求mu lambda  
Mu = [0 for x in range(n-2)]  
Lambda = [0 for x in range(n-2)]  
h = [0 for x in range(n-1)]  
for i in range(n-1):  
    h[i] = x_n[i+1] - x_n[i]  
for i in range(n-2):  
    Mu[i] = h[i] / (h[i] + h[i+1])  
    Lambda[i] = 1 - Mu[i]  
  
# 追赶法求解线性方程组  
c = [1] + Lambda  
a = Mu + [1]  
b = 2  
# 消元  
beta = [b]  
y = [d[0]]  
for i in range(1, n):  
    l_ = a[i-1] / beta[i-1]  
    beta.append(2 - l_ * c[i-1])  
    y.append(d[i] - l_ * y[i-1])  
# 回代  
M = [y[n-1]/beta[n-1]]  
for i in range(n-2, -1, -1):  
    M.insert(0, (y[i] - c[i] * M[0])/beta[i])  
  
# 输出S(x)  
Sx = []  
for i in range(0, n-1):  
    x = symbols('x')  
    SX = y_n[i] + (f1[i] - (1/3 * M[i] + 1/6 * M[i+1]) * h[i]) * (x - x_n[i]) + \  
        1/2 * M[i] * ((x - x_n[i]) ** 2) + 1/(6 * h[i]) * (M[i+1] - M[i]) * ((x - x_n[i]) ** 3)  
    SX = expand(SX)  
    Sx.append(SX)  
    print(SX, end="")  
    print("   ({:d}<x<{:d})".format(int(x_n[i]), int(x_n[i+1])))  
  
  
# 插值  
def S(value):  
    for j in range(0, n-1):  
        if (value >= x_n[j]) and (value < x_n[j+1]):  
            return Sx[j].subs(x, value)  
    return "不在插值范围内"  
  
  
print("i     x     S(i+0.5)")  
for k in range(0, n-1):  
    dot = float(k + 0.5)  
    print("{:d}    {:.1f}    {:.5f}".format(k, dot, S(dot)))

4.3 程序运行结果

4.3.1 程序输入

请输入插值点个数：11
请输入插值节点的横坐标，以空格分隔：0 1 2 3 4 5 6 7 8 9 10
请输入插值节点的函数值，以空格分隔：2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80
请输入左边界条件:y'(x0)=0.8
请输入右边界条件:y'(xn)=0.2

• 通过键盘输入的方式将第一型3次样条插值函数的所需参数进行输入，增强程序的泛用性

4.3.2 3次样条插值函数

-0.00851403685003562*x**3 - 0.00148596314996439*x**2 + 0.8*x + 2.51   (0<x<1)
-0.00445788944989296*x**3 - 0.0136544053503924*x**2 + 0.812168442200428*x + 2.50594385259986   (1<x<2)
-0.0036544053503928*x**3 - 0.0184753099473933*x**2 + 0.82181025139443*x + 2.49951597980386   (2<x<3)
-0.0409244891485364*x**3 + 0.316955444235899*x**2 - 0.184482011155446*x + 3.50580824235373   (3<x<4)
0.10735236194454*x**3 - 1.46236676888101*x**2 + 6.9328068413122*x - 5.98391022760313   (4<x<5)
-0.268484958629623*x**3 + 4.17519303973142*x**2 - 21.25499220175*x + 40.9957548441672   (5<x<6)
0.426587472573951*x**3 - 8.33611072193291*x**2 + 53.812830368236*x - 109.139890295805   (6<x<7)
-0.267864931666182*x**3 + 6.24738976710989*x**2 - 48.2716730550636*x + 129.057284358561   (7<x<8)
0.0548722540907767*x**3 - 1.49830269105712*x**2 + 13.6938666102725*x - 36.184154749002   (8<x<9)
0.0583759153030747*x**3 - 1.59290154378917*x**2 + 14.5452562848609*x - 38.7383237727672   (9<x<10)

整理得
$$\begin{array}{rl}S(x)& =-0.00851x^3-0.00149x^2+0.8x+2.51(0\le x<1)\\ S(x)& =-0.00446x^3-0.01365x^2+0.81217x+2.50594(1\le x<2)\\ S(x)& =-0.00365x^3-0.01848x^2+0.82181x+2.49952(2\le x<3)\\ S(x)& =-0.04092x^3+0.31696x^2-0.18448x+3.50581(3\le x<4)\\ S(x)& =0.10735x^3-1.46237x^2+6.93281x-5.98391(4\le x<5)\\ S(x)& =-0.26848x^3+4.17519x^2-21.25499x+40.99575(5\le x<6)\\ S(x)& =0.42659x^3-8.33611x^2+53.81283x-109.13989(6\le x<7)\\ S(x)& =-0.26786x^3+6.24739x^2-48.27167x+129.05728(7\le x<8)\\ S(x)& =0.05487x^3-1.49830x^2+13.69387x-36.18415(8\le x<9)\\ S(x)& =0.05838x^3-1.59290x^2+14.54526x-38.73832(9\le x<10)\end{array}$$

4.3.3 非插值点插值结果

i     x     S(i+0.5)
0    0.5    2.90856
1    1.5    3.67843
2    2.5    4.38147
3    3.5    4.98819
4    4.5    5.38328
5    5.5    5.72370
6    6.5    5.59441
7    7.5    5.42989
8    8.5    5.65977
9    9.5    5.73230

4.4 总结感悟

• 3次样条插值函数的计算编程实现相当的复杂，其中还涉及到线性方程组求解的编程，但是由于第一型3次样条插值函数的三弯矩方程可以使用追赶法进行编程求解，省去了Guass消元法编程的麻烦。
• 该Python编程主要采用的是List编程实现，其中存在一定的冗余编程，如过多的使用了for循环，如果考虑采用numpy数组编程的话，引入一些数据处理与矩阵运算的函数可在一定程度上简化编程，但是违背了上机题一步步底层实现的初衷。
• 3次样条插值函数是一种广泛应用于工程实践当中的函数插值方法，具有一些优良的性质，如2阶连续可导，该问题即为汽车门曲线的插值求解，对我们后续的工程实践具有很强的指导意义。