数值分析实验作业（一）——三次样条插值

前言

数值分析课的第一次实验作业。在此浅作记录。

问题描述

编程实现三次样条插值方法。

问题分析

按照《数值分析》课本第43~44页~~（懒得找了）~~建立样条插值函数的方法，先输入节点的数量和边界条件的种类，再依次输入各节点的值及其对应的函数值。根据这些条件确定$h_j$、$\mu_j$、$\lambda_j$与$d_j$的值，再运用矩阵求解得到各$M_j$的值，最后便可计算出各区间的$S(x)$的表达式。

在进行数据存储与矩阵运算时引入了Python中的numpy库并使用其相应功能，这样就不用手写矩阵运算了。

完整代码

import numpy as np
import math
from numpy.linalg import (solve)

m = int(input("输入节点的数量: "))
n = m-1                            #使后续程序的下标与课本内保持一致
X = np.zeros((m))                #节点值
Y = np.zeros((m))                #对应的函数值
H = np.zeros((m-1))
Mu = np.zeros((m))
Lambda = np.zeros((m))
D = np.zeros((m, 1))
A = np.zeros((m, m))
S = np.zeros((n, 4))            #各区间对应表达式的k（k=0,1,2,3）次项的系数


def Deviation(i, j):                        #计算f([x_i, x_j])的均差
    return (Y[j]-Y[i])/(X[j]-X[i])


opt = input("输入边界条件种类(1~3): ")

for i in range(n+1):
    x, y = input("输入节点x_%d的值及对应的函数值: " % i).split()
    X[i] = float(x)
    Y[i] = float(y)
for j in range(n):
    H[j] = X[j+1]-X[j]

if opt == "1" or opt == "2":
    diff_0, diff_n = input("输入两个边界点的"+opt+"阶导的值: ").split()
    diff_0 = float(diff_0)
    diff_n = float(diff_n)
    for j in range(1, n):
        Mu[j] = H[j-1]/(H[j-1]+H[j])
        Lambda[j] = 1.0-Mu[j]
        D[j] = 6*(Deviation(j, j+1)-Deviation(j-1, j))/(H[j-1]+H[j])
    if opt == "1":                            #第一种边界条件，即已知两端的一阶导数值
        Lambda[0] = Mu[n] = 1
        D[0] = 6*(Deviation(0, 1)-diff_0)/H[0]
        D[n] = 6*(diff_n-Deviation(n-1, n))/H[n-1]
    else:                                    #第二种边界条件，即已知两端的二阶导数值
        Lambda[0] = Mu[n] = 0
        D[0] = 2*diff_0
        D[n] = 2*diff_n
    for j in range(n):
        A[j][j] = 2
        A[j][j+1] = Lambda[j]
        A[j+1][j] = Mu[j+1]
    A[n][n] = 2
    M = solve(A, D)                            #求解矩阵表达式AM = D的解
    for j in range(n):
        S[j][3] = (M[j+1]-M[j])/(6*H[j])
        S[j][2] = (M[j]*X[j+1]-M[j+1]*X[j])/(2*H[j])
        S[j][1] = (M[j+1]*math.pow(X[j], 2)-M[j]*math.pow(X[j+1], 2) +
                   2*(Y[j+1]-Y[j]))/(2*H[j])+(M[j]-M[j+1])*H[j]/6
        S[j][0] = (M[j]*math.pow(X[j+1], 3)-M[j+1]*math.pow(X[j], 3)+6 *
                   (Y[j]*X[j+1]-Y[j+1]*X[j]))/(6*H[j])+(M[j+1]*X[j]-M[j]*X[j+1])*H[j]/6
    print("S(x)")
    for j in range(n):
        print("= (%f*x^3) + (%f*x^2) + (%f*x) + (%f), %f <= x <= %f" %
              (S[j][3], S[j][2], S[j][1], S[j][0], X[j], X[j+1]))

elif opt == "3":                            #第三种边界条件，即已知该函数为周期函数
    D_temp = np.zeros((m-1, 1))                #用于矩阵求解
    A_temp = np.zeros((m-1, m-1))            #用于矩阵求解
    for j in range(1, n):
        Mu[j] = H[j-1]/(H[j-1]+H[j])
        Lambda[j] = 1.0-Mu[j]
        D[j] = 6*(Deviation(j, j+1)-Deviation(j-1, j))/(H[j-1]+H[j])
    Lambda[n] = H[0]/(H[n-1]+H[0])
    Mu[n] = 1.0-Lambda[n]
    D[n] = 6*(Deviation(0, 1)-Deviation(n-1, n))/(H[n-1]+H[0])
    for j in range(n):
        D_temp[j] = D[j+1]
    for j in range(n-1):
        A_temp[j][j] = 2
        A_temp[j][j+1] = Lambda[j+1]
        A_temp[j+1][j] = Mu[j+2]
    A_temp[n-1][n-1] = 2
    A_temp[0][n-1] = Mu[1]
    A_temp[n-1][0] = Lambda[n]
    M_temp = solve(A_temp, D_temp)
    M = np.zeros((m, 1))
    for j in range(n):
        M[j+1] = M_temp[j]                    #保持下标与课本内一致
    M[0] = M[n]
    for j in range(n):
        S[j][3] = (M[j+1]-M[j])/(6*H[j])
        S[j][2] = (M[j]*X[j+1]-M[j+1]*X[j])/(2*H[j])
        S[j][1] = (M[j+1]*math.pow(X[j], 2)-M[j]*math.pow(X[j+1], 2) +
                   2*(Y[j+1]-Y[j]))/(2*H[j])+(M[j]-M[j+1])*H[j]/6
        S[j][0] = (M[j]*math.pow(X[j+1], 3)-M[j+1]*math.pow(X[j], 3)+6 *
                   (Y[j]*X[j+1]-Y[j+1]*X[j]))/(6*H[j])+(M[j+1]*X[j]-M[j]*X[j+1])*H[j]/6
    difference = float(X[n]-X[0])            #该函数的周期
    print("S(x)")
    for j in range(n):
        print("= (%f*x^3) + (%f*x^2) + (%f*x) + (%f), %f+k*%f <= x <= %f+k*%f, k为整数" %
              (S[j][3], S[j][2], S[j][1], S[j][0], X[j], difference, X[j+1], difference))

后记

毕竟已经有现成的公式可以直接抄了，所以写起来思路还是比较清晰的。就是变量太多，写着写着自己还搞混过几次，怎么会事呢.jpg