一、Lagrange插值
L n ( x ) = ∑ i = 0 n l i ( x ) f ( x i ) L_n(x)=\sum_{i=0}^{n}l_i(x)f(x_i) Ln(x)=∑i=0nli(x)f(xi),其中 l i ( x ) = ∏ j = 0 , j ≠ i n x − x j x i − x j l_i(x)=\prod_{j=0,j\ne i}^n \frac{x-x_j}{x_i-x_j} li(x)=∏j=0,j=inxi−xjx−xj。
# -*- 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 ) + ⋯ + f [ x 0 , ⋯ , x n ] ( x − x 0 ) ⋯ ( x − x n − 1 ) 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}) Nn(x)=f[x0]+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,⋯,xn](x−x0)⋯(x−xn−1),其中 f [ x 0 , ⋯ , x k ] = f [ x 1 , ⋯ , x k ] − f [ x 0 , ⋯ , x k − 1 ] x k − x 0 , k = 1 , ⋯ , n 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 f[x0,⋯,xk]=xk−x0f[x1,⋯,xk]−f[x0,⋯,xk−1],k=1,⋯,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
n表示函数次幂,代码中的
n
n
n表示点数,一般次幂=点数-1。
三、两点三次Hermite插值
H 3 ( x ) = ( 1 + 2 x − x k x k + 1 − x k ) ( x − x k + 1 x k − x k + 1 ) 2 y k + ( 1 + 2 x − x k + 1 x k − x k + 1 ) ( x − x k x k + 1 − x k ) 2 y k + 1 + ( x − x k ) ( x − x k + 1 x k − x k + 1 ) 2 m k + ( x − x k + 1 ) ( x − x k x k + 1 − x k ) 2 m k + 1 H_3(x)=(1+2\frac{x-x_k}{x_{k+1}-x_{k}})(\frac{x-x_{k+1}}{x_k-x_{k+1}})^2y_k+(1+2\frac{x-x_{k+1}}{x_{k}-x_{k+1}})(\frac{x-x_{k}}{x_{k+1}-x_{k}})^2y_{k+1}+(x-x_k)(\frac{x-x_{k+1}}{x_k-x_{k+1}})^2m_k+(x-x_{k+1})(\frac{x-x_k}{x_{k+1}-x_k})^2m_{k+1} H3(x)=(1+2xk+1−xkx−xk)(xk−xk+1x−xk+1)2yk+(1+2xk−xk+1x−xk+1)(xk+1−xkx−xk)2yk+1+(x−xk)(xk−xk+1x−xk+1)2mk+(x−xk+1)(xk+1−xkx−xk)2mk+1,其中 m k = f ′ ( x k ) m_k=f'(x_k) mk=f′(xk)
# -*- 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_)
四、三次样条插值
μ
j
M
j
−
1
+
2
M
j
+
λ
j
M
j
+
1
=
d
j
,
j
=
1
,
2
,
⋯
,
n
−
1
\mu_jM_{j-1}+2M_j+\lambda_jM_{j+1}=d_j,j=1,2,\cdots,n-1
μjMj−1+2Mj+λjMj+1=dj,j=1,2,⋯,n−1,其中
μ
j
=
h
j
−
1
h
j
−
1
+
h
j
,
λ
j
=
h
j
h
j
−
1
+
h
j
,
d
j
=
6
f
[
x
j
,
x
j
+
1
]
−
f
[
x
j
−
1
,
x
j
]
h
j
−
1
+
h
j
=
6
f
[
x
j
−
1
,
x
j
,
x
j
+
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}]
μj=hj−1+hjhj−1,λj=hj−1+hjhj,dj=6hj−1+hjf[xj,xj+1]−f[xj−1,xj]=6f[xj−1,xj,xj+1];
对于第一类边界条件,还有两个方程:
2
M
0
+
M
1
=
6
h
0
(
f
[
x
0
,
x
1
]
−
f
0
′
)
;
M
n
−
1
+
2
M
n
=
6
h
n
−
1
(
f
n
′
−
f
[
x
n
−
1
,
x
n
]
)
2M_0+M_1=\frac{6}{h_0}(f[x_0,x_1]-f'_0);M_{n-1}+2M_n=\frac{6}{h_{n-1}}(f'_n-f[x_{n-1},x_n])
2M0+M1=h06(f[x0,x1]−f0′);Mn−1+2Mn=hn−16(fn′−f[xn−1,xn])。
S
(
x
)
=
M
j
(
x
j
+
1
−
x
)
3
6
h
j
+
M
j
+
1
(
x
−
x
j
)
3
6
h
j
+
(
y
j
−
M
j
h
j
2
6
)
x
j
+
1
−
x
h
j
+
(
y
j
+
1
−
M
j
+
1
h
j
2
6
)
x
−
x
j
h
j
,
j
=
0
,
1
,
⋯
,
n
−
1
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_jh_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
S(x)=Mj6hj(xj+1−x)3+Mj+16hj(x−xj)3+(yj−6Mjhj2)hjxj+1−x+(yj+1−6Mj+1hj2)hjx−xj,j=0,1,⋯,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_)