联系数计加油站获取论文和代码等资源。
邮箱:MathComputerSharing@outlook.com
微信公众号:数计加油站
致谢:清风数学建模
什么是插值😕
插值，顾名思义，就是插进几个值。往哪里插 ？往已知的数据中间插。
无论是做实验，还是做记录，我们很难获取到连续的数据值。我们获取到的数据，往往都是一个个离散的数据点。
相邻的两个数据点之间的数据，到底是怎么样的呢?我们不能断定。但可以通过插值去“无中生有”地去估计。
插值的较为正式的定义是:寻找到一个函数，它穿过所有已知数据点。 通俗的讲，就是要找到一根线，去串起几个已知的点。
多项式函数插值与龙格现象
多项式函数可以说是最简单最常见的函数了吧，谁赞成，谁反对❓
多项式函数不仅性质简单，而且易于在计算机系统中表示和运算，所以常常是人们迫害的对象。
可以证明， n + 1 n+1 n+1个数据点可以唯一确定一个 n n n元多项式函数。
但这样直白的插值，容易出现问题。 n n n太大的话，由于高幂方运算，有时候多项式函数的图像会上下大幅度震荡，x微小的改变可能导致y的大幅度改变。这种现象称为“龙格震荡现象”。虽然我们并不知道原始数据点满足的函数是不是也这么洒脱😇 ，但是一般来说，出现了龙格现象的多项式函数插值的效果是不太好的。
分段插值:越细越准
解决单个多项式函数容易出现龙格现象的方法，其实很简单。那就是每几个点用次数低的多项式函数插一下，最后合成一个整的分段函数。
比如最常见也是最简单的分段插值：分段线性插值——把每两个点用直线连一下。
事实上，很多函数绘图软件也是用分段线性插值的思想实现的。你给它一串数据点，它连一下，只要数据点的个数够多，最后的函数图像就越弯越像一条名副其实的函数曲线。也就是说，分段线性插值，数据点越多，就越接近于真实规律。
样条插值:弯弯的更好看💙
分段线性插值在连接点太尖了。我们都知道，太尖的函数不太好，它在尖的地方不存在导数。
k k k次样条插值要求插出来的整体函数有 k − 1 k-1 k−1阶连续导数。常用的是三次样条函数。
样条插值插出来的函数特别得光滑，特别的好看。而我们遇到的函数规律，一般都是弯弯的。
一篇简单的小论文
Matlab代码
待定系数方程组法
function y = SolveInterpolate(x0,y0,x)
% 使用方程组法得到单多项式插值函数
% x0,y0:已知的原始数据点
% x,y:插值得到的点
% x,y,x0,y0均为列向量
%% 输入参数解方程 X a = y0
d = vander(x0);
a = d \ y0;
y = polyval(a,x);
end
拉格朗日插值法
function y = LagrangeInterpolate(x0,y0,x)
% 使用拉格朗日插值法得到单多项式插值函数
% x0,y0:已知的原始数据点
% x,y:插值得到的点
%% 拉格朗日插值
n = length(x0);
m = length(x);
y = ones(m,1);
for k = 1:m
x_new = x(k);
y_new = 0;
% 构造第i个基函数l_i,并代入第k个新x值,对l_i进行累加得到y_new
for i = 1:n
l_i=1;
for j = 1:n
if j~=i
l_i = l_i * (x_new-x0(j))/(x0(i)-x0(j));
end
end
y_new = y_new + y0(i)*l_i;
end
y(k) = y_new;
end
end
求差商表
function dqTable = DifferQuotient(x,y)
% 根据输入的x,y(n+1维)计算n阶差商
n = length(x)-1;
dqTable = cell(n+1,1);
dqTable{1} = y;
for j = 1:n
delta_y = diff(dqTable{j});
for i = 1:n-j+1
delta_x = x(i+j)-x(i);
end
dqTable{j+1} = delta_y ./ delta_x;
end
dqTable = dqTable(2:end,:);
end
牛顿插值
function y = NewtonInterpolate(x0,y0,x)
% 使用牛顿插值法得到单多项式插值函数
% x0,y0:已知的原始数据点
% x,y:插值得到的点
%% 牛顿插值
n = length(x0);
m = length(x);
y = ones(m,1);
% 计算差商
dqTable = DifferQuotient(x0,y0);
% 求第i个基函数并代入第k个新x值，结果相加
for k=1:m
x_new = x(k);
Product = 1;
N = y0(1);
for i=1:n-1
Product = Product*(x_new-x0(i));
N = N + dqTable{i}(1,1)*Product;
end
y(k,1) = N;
end
end
插值问题主体
clear,clc
%% 定义原函数
f = @(x) 1./(1+x.^2);
%% 绘制原函数
x = -5:0.01:5;
y = f(x);
plot(x,y,'black');
hold on;
%% 绘制单段多项式插值函数
color = ["green","blue","red"];
for n = [6,8,10]
x0 = linspace(-5,5,n+1);
y0 = f(x0);
new_y = NewtonInterpolate(x0,y0,x);
plot(x0,y0,color(n/2-2)+'o',x,new_y,color(n/2-2));
hold on
end
legend("原函数","","n=6","","n=8","","n=10","Location","south");
clc,clear
%% 定义原函数
f = @(x) 1./(1+x.^2);
%% 绘制原函数
x = -5:0.01:5;
y = f(x);
subplot(221)
plot(x,y,'black');
xlabel("原函数")
%% 绘制线性插值函数
color = ["green","blue","red"];
for n = [8,16,24]
x0 = linspace(-5,5,n+1);
y0 = f(x0);
new_y = interp1(x0,y0,x,"linear");
subplot(221+n/8);
plot(x0,y0,color(n/8)+'x',x,new_y,color(n/8));
xlabel("n="+n);
end
clc,clear
%% 定义原函数
f = @(x) 1./(1+x.^2);
%% 绘制原函数
x = -5:0.01:5;
y = f(x);
subplot(221)
plot(x,y,'black');
xlabel("原函数")
%% 绘制三次样条插值函数
color = ["green","blue","red"];
for n = [6,12,18]
x0 = linspace(-5,5,n+1);
y0 = f(x0);
new_y = interp1(x0,y0,x,"spline");
subplot(221+n/6);
plot(x0,y0,color(n/6)+'x',x,new_y,color(n/6));
xlabel("n="+n);
end
Python代码
初始化
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
np.set_printoptions(suppress=True)
plt.rcParams["font.sans-serif"] = ["SimHei"] #设置字体
plt.rcParams["axes.unicode_minus"] = False # 解决图像中的“-”负号的乱码问题
待定系数方程法
def solve_interpolate(x0, y0, x):
"""使用方程组法得到单段多项式插值函数"""
# 输入参数解方程 X a = y0
d = np.vander(x0)
a = np.dot(np.linalg.inv(d), y0)
y = np.polyval(a, x)
return y
拉格朗日插值法
def langrange_interpolate(x0, y0, x):
"""使用拉格朗日插值法得到单多项式插值函数"""
# 参数初始化
n = len(x0)
m = len(x)
y = np.ones(m)
for k in range(m):
x_new = x[k]
y_new = 0
# 构造第i个基函数l_i,并代入第k个新x值,对l_i进行累加得到y_new
for i in range(n):
l_i = 1
for j in range(n):
if j != i:
l_i = l_i * (x_new - x0[j]) / (x0[i] - x0[j])
y_new = y_new + y0[i] * l_i
y[k] = y_new
return y
牛顿插值法
def differ_quotient(x, y):
"""根据输入的x,y(n+1维)计算n阶差商"""
# 参数初始化
n = len(x) - 1
dqTable = [y]
# 求差商表
for j in range(1, n + 1):
delta_y = np.diff(dqTable[j - 1])
delta_x = 0
for i in range(1, n - j + 2):
delta_x = x[i + j - 1] - x[i - 1]
dqTable.append(np.array(delta_y / delta_x))
dqTable = dqTable[1:]
return dqTable
def newton_interpolate(x0, y0, x):
"""使用牛顿插值法得到单多项式插值函数"""
# 参数初始化
n = len(x0)
m = len(x)
y = np.ones(m)
dqTable = differ_quotient(x0, y0)
# 求第i个基函数并代入第k个新x值，结果相加
for k in range(m):
x_new = x[k]
Product = 1
N = y0[0]
for i in range(n - 1):
Product = Product * (x_new - x0[i])
N = N + dqTable[i][0] * Product
y[k] = N
return y
龙格震荡现象
# 定义原函数
f = lambda x: 1 / (1 + x ** 2)
# 绘制原函数
x = np.linspace(-5, 5, 1000)
y = f(x)
fig, ax = plt.subplots()
ax.plot(x, y, color="black")
# 绘制单段多项式插值函数
color = ["green", "blue", "red"]
labels = ["n=6", "n=8", "n=10"]
line_width = 0.5
for n in [6, 8, 10]:
x0 = np.linspace(-5, 5, n + 1)
y0 = f(x0)
new_y = newton_interpolate(x0, y0, x)
ax.scatter(x0, y0, color=color[int(n / 2 - 2) - 1], marker="o", facecolor="white", linewidth=line_width)
ax.plot(x, new_y, color=color[int(n / 2 - 2) - 1], linewidth=line_width, label=labels[int(n / 2 - 2) - 1])
plt.legend()
分段线性插值
# 定义原函数
f = lambda x: 1 / (1 + x ** 2)
# 绘制原函数
x = np.linspace(-5, 5, 1000)
y = f(x)
plt.subplot(2,2,1)
plt.plot(x, y, color="black")
plt.xlabel("原函数")
# 绘制分段线性插值函数
color = ["green", "blue", "red"]
line_width = 0.5
for n in [8, 16, 24]:
x0 = np.linspace(-5, 5, n + 1)
y0 = f(x0)
linear_interp = interpolate.interp1d(x0,y0,'linear')
new_y = linear_interp(x)
plt.subplot(2,2,int(n/8+1))
plt.scatter(x0, y0, color=color[int(n / 8) - 1], marker="o", facecolor="white", linewidth=line_width)
plt.plot(x, new_y, color=color[int(n /8) - 1], linewidth=line_width)
plt.xlabel(f"$n={n}$")
plt.subplots_adjust(wspace=0.3,hspace=0.3)
三次样条插值
# 定义原函数
f = lambda x: 1 / (1 + x ** 2)
# 绘制原函数
x = np.linspace(-5, 5, 1000)
y = f(x)
plt.subplot(2,2,1)
plt.plot(x, y, color="black")
plt.xlabel("原函数")
# 绘制三次样条插值函数
color = ["green", "blue", "red"]
line_width = 0.5
for n in [6, 12, 18]:
x0 = np.linspace(-5, 5, n + 1)
y0 = f(x0)
spline_interp = interpolate.CubicSpline(x0,y0)
new_y = spline_interp(x)
plt.subplot(2,2,int(n/6+1))
plt.scatter(x0, y0, color=color[int(n / 6) - 1], marker="o", facecolor="white", linewidth=line_width)
plt.plot(x, new_y, color=color[int(n /6) - 1], linewidth=line_width)
plt.xlabel(f"$n={n}$")
plt.subplots_adjust(wspace=0.3,hspace=0.3)
7741
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
