插值篇(数值分析)

一、Lagrange插值

$L_n(x)={\sum }_{i=0}^n{l}_i(x)f(x_i)$，其中$l_i(x)={\prod }_{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}$。

# -*- coding: utf-8 -*-
"""
Created on Tue Mar 15 09:57:38 2022

@author: 86155
"""
import numpy as np
import matplotlib.pyplot as plt

# 计算基函数
def l_i(i, x, x0):
    den = 1 # 分母
    mol = 1 # 分子
    for j in range(len(x)):
        if j != i:
            den *= x[i] - x[j]
            mol *= x0 - x[j]
    return mol/den
# 计算插值结果
def L_n(n, x, y, x0):
    result = 0
    for i in range(n):
        result += l_i(i, x, x0)*y[i]
    return result

# 插值点
x = np.array([-2,0,1,2])
y = np.array([17, 1, 2, 17])
plt.scatter(x, y)
# 计算插值函数
a = np.min(x)
b = np.max(x)
n = len(x)
x_ = np.linspace(a, b, 20)
y_ = [] # y的近似值
for i in range(len(x_)):
    y_.append(L_n(n, x, y, x_[i]))
y_ = np.array(y_)
plt.plot(x_, y_)

二、Newton插值

$N_n(x)=f[x_0]+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdots +f[x_0,\cdots ,x_n](x-x_0)\cdots (x-x_{n-1})$，其中$f[x_0,\cdots ,x_k]=\frac{f[x_1,\cdots ,x_k]-f[x_0,\cdots ,x_{k-1}]}{x_k-x_0},k=1,\cdots ,n$

# -*- coding: utf-8 -*-
"""
Created on Wed Mar 16 21:22:18 2022

@author: 86155
"""
import numpy as np
import matplotlib.pyplot as plt


# 计算差商
def f(x, y, i, j):
    if j - i > 1:
        return (f(x, y, i+1, j) - f(x, y, i, j-1))/(x[j] - x[i])
    else:
        return (y[j] - y[i])/(x[j] - x[i])

# 计算结果
def N_n(x, y, x0):
    n = len(x)
    result = f(x, y, 0, n-1)*(x0 - x[n-2])
    for i in range(n-2, 0, -1):
        result += f(x, y, 0, i)
        result *= (x0 - x[i-1])
        print(i)
        print(result)
    result += y[0]
    return result


# 插值点
x = np.array([0.4,0.55,0.65,0.80, 0.90])
y = np.array([0.41075, 0.57815, 0.69675, 0.88811, 1.02652])
plt.scatter(x, y)
# 计算插值函数
a = np.min(x)
b = np.max(x)
n = len(x)
x_ = np.linspace(a, b, 20)
y_ = [] # y的近似值
for i in range(len(x_)):
    y_.append(N_n(x, y, x_[i]))
y_ = np.array(y_)
plt.plot(x_, y_)

注意：公式中的$n$表示函数次幂，代码中的$n$表示点数，一般次幂=点数-1。

三、两点三次Hermite插值

$H_3(x)=(1+2\frac{x-x_k}{x_{k+1}-x_k})(\frac{x-x_{k+1}}{x_k-x_{k+1}}{)}^2{y}_k+(1+2\frac{x-x_{k+1}}{x_k-x_{k+1}})(\frac{x-x_k}{x_{k+1}-x_k}{)}^2{y}_{k+1}+(x-x_k)(\frac{x-x_{k+1}}{x_k-x_{k+1}}{)}^2{m}_k+(x-x_{k+1})(\frac{x-x_k}{x_{k+1}-x_k}{)}^2{m}_{k+1}$，其中$m_k=f^{\prime }(x_k)$

# -*- coding: utf-8 -*-
"""
Created on Wed Mar 23 22:52:06 2022

@author: 86155
"""
## 两点三次Hermite插值
import numpy as np
import matplotlib.pyplot as plt

def seg_fun(x0, x, y, dy):
    f1 = (1+2*(x0-x[0])/(x[1]-x[0]))*((x0-x[1])/(x[0]-x[1]))**2*y[0]
    f2 = (1+2*(x0-x[1])/(x[0]-x[1]))*((x0-x[0])/(x[1]-x[0]))**2*y[1]
    f3 = (x0-x[0])*((x0-x[1])/(x[0]-x[1]))**2*dy[0]
    f4 = (x0-x[1])*((x0-x[0])/(x[1]-x[0]))**2*dy[1]
    return f1+f2+f3+f4

x = [1/4, 1, 9/4]
y = [1/8, 1, 27/8]
dy = [3/4, 3/2, 9/4]
N = 100
x_ = []
y_ = []
for i in range(len(x)-1):
    h = (x[i+1] - x[i])/N
    for j in range(N):
        x0 = x[i] + j*h
        x_.append(x0)
        y_.append(seg_fun(x0, [x[i], x[i+1]], [y[i], y[i+1]], [dy[i], dy[i+1]]))
x = np.array(x)
y = np.array(y)
plt.scatter(x, y)
x_ = np.array(x_)
y_ = np.array(y_)
plt.plot(x_, y_)

四、三次样条插值

${\mu }_j{M}_{j-1}+2M_j+{\lambda }_j{M}_{j+1}=d_j,j=1,2,\cdots ,n-1$，其中${\mu }_j=\frac{h_{j-1}}{h_{j-1}+h_j},{\lambda }_j=\frac{h_j}{h_{j-1}+h_j},d_j=6\frac{f[x_j,x_{j+1}]-f[x_{j-1},x_j]}{h_{j-1}+h_j}=6f[x_{j-1},x_j,x_{j+1}]$；
对于第一类边界条件，还有两个方程：$2M_0+M_1=\frac{6}{h_0}(f[x_0,x_1]-f_0^{\prime });M_{n-1}+2M_n=\frac{6}{h_{n-1}}(f_n^{\prime }-f[x_{n-1},x_n])$。
$S(x)=M_j\frac{(x_{j+1}-x{)}^3}{6h_j}+M_{j+1}\frac{(x-x_j{)}^3}{6h_j}+\left(y_j-\frac{M_j{h}_j^2}{6}\right)\frac{x_{j+1}-x}{h_j}+\left(y_{j+1}-\frac{M_{j+1}{h}_j^2}{6}\right)\frac{x-x_j}{h_j},j=0,1,\cdots ,n-1$

# -*- coding: utf-8 -*-
"""
Created on Tue Mar 29 10:03:39 2022

@author: 86155
"""
import numpy as np
import matplotlib.pyplot as plt


# 计算差商
def f(x, y, i, j):
    if j - i > 1:
        return (f(x, y, i+1, j) - f(x, y, i, j-1))/(x[j] - x[i])
    else:
        return (y[j] - y[i])/(x[j] - x[i])

# 分段函数
def seg_fun(x, y, M, i, x0):
    f1 = M[i][0]*(x[i+1] - x0)**3/(6*h[i])
    f2 = M[i+1][0]*(x0 - x[i])**3/(6*h[i])
    f3 = (y[i] - M[i][0]*h[i]**2/6)*(x[i+1]-x0)/h[i]
    f4 = (y[i+1] - M[i+1][0]*h[i]**2/6)*(x0 - x[i])/h[i]
    return f1+f2+f3+f4

x = [27.7, 28, 29, 30]
y = [4.1, 4.3, 4.1, 3.0]
f0_ = 3.0
fn_ = -4.0 # 第一类边界条件
n = len(x)
N = 1000 # 一个间隔划分的区间段

h = [x[i+1] - x[i] for i in range(n - 1)]
mu = [h[i-1]/(h[i-1]+h[i]) for i in range(1, len(h))] # 下次对角线
mu.append(1)
con = [2 for i in range(n)] # 主对角线
lam = [h[i]/(h[i-1]+h[i]) for i in range(1, len(h))] # 上次对角线
lam.insert(0, 1)
A = np.zeros((n, n))
d = np.zeros((n, 1))
for i in range(n):
    if i == 0:
        A[i][i] = con[i]
        A[i][i+1] = lam[i]
        d[i][0] = 6/h[0]*(f(x, y, 0, 1) - f0_)
    elif i == n-1:
        A[i][i-1] = mu[i-1]
        A[i][i] = con[i]
        d[i][0] = 6/h[n-2]*(fn_ - f(x, y, n-2, n-1))
    else:
        A[i][i-1] = mu[i-1]
        A[i][i] = con[i]
        A[i][i+1] = lam[i]
        d[i][0] = 6*f(x, y, i-1, i+1)

M = np.dot(np.linalg.inv(A),d)

x_ = []
y_ = []
for i in range(n - 1):
    h0 = h[i]/N
    for j in range(N):
        x0 = x[i] + j*h0
        x_.append(x0)
        y_.append(seg_fun(x, y, M, i, x0))
x = np.array(x)
y = np.array(y)
plt.scatter(x, y)
x_ = np.array(x_)
y_ = np.array(y_)
plt.plot(x_, y_)