参考书籍:数值分析 第五版 李庆杨 王能超 易大义编 第2章 插值法
文章声明:如有发现错误,欢迎批评指正
2.1引言
2.1.1插值问题的提出
一些概念:插值函数,插值节点,插值区间,插值法,插值多项式,多项式插值,分段插值,三角插值。
2.1.2多项式插值
范德蒙矩阵的性质,定理1:存在唯一。
补充1:多项式插值
题面描述:任意给出n个节点,求出多项式插值函数并且绘制多项式插值函数曲线。
Python
注释:不能支持绘制多条函数曲线
import matplotlib.pyplot as plt
import numpy as np
def this_function(x):
cnt=0
for i in range(len(A)):
cnt+=A[i]*pow(x,i)
return cnt
print("注意输入的x与y一一对应。")
ltx=[eval(i) for i in input("输入n个节点的x,使用空格分割:").strip().split()]
lty=[eval(i) for i in input("输入n个节点的y,使用空格分割:").strip().split()]
lt=[]
for i in range(len(ltx)):
lt.append([pow(ltx[i],j) for j in range(len(ltx))])
lt=np.array(lt);lty=np.array(lty)
A=np.linalg.solve(lt,lty)
x=[0+i*0.001 for i in range(10001)]
y=[this_function(i) for i in x]
plt.figure(figsize=(16,9))
plt.scatter(ltx,lty,s=10,c="r")
plt.plot(x,y,linewidth=1,color="b")
plt.grid();plt.show()
Matlab
function []=duoxiangshichazhi()
function [cnt]=this_function(x)
cnt=0;
for i=1:length(A)
cnt=cnt+A(i)*power(x,i-1);
end
end
disp('注意输入的x与y一一对应。');
ltx=input('输入n个节点的x,使用空格分割:','s');
ltx=str2num(ltx);
lty=input('输入n个节点的y,使用空格分割:','s');
lty=str2num(lty);
lt=zeros(length(ltx));
for i=1:length(ltx)
lt(i,:)=power(ltx(i),[0:length(ltx)-1]);
end
A=(inv(lt)*(lty'));
x=[0:0.001:10];
y=arrayfun(@this_function,x);
figure('Position',[0 0 1600 1200]);
hold on;
scatter(ltx,lty,10,'r','filled');
plot(x,y,'b','LineWidth', 1);
grid on;
hold off;
end
结果
输入
输出
注意事项:这里必须求得矩阵的逆但是博主不会只好掉包(求逆的包)。接下来的任何方法不会掉包,我们从底层来构建。
2.2拉格朗日插值
2.2.1线性插值与抛物线插值
n个节点:n-1次插值基函数
2.2.2拉格朗日插值多项式
定义1:n次插值基函数,拉格朗日插值多项式
2.2.3插值余项与误差估计
插值多项式的余项,定理2:插值余项具体是个什么
补充2:拉格朗日插值
题面描述:任意给出n个节点,求出拉格朗日插值函数并且绘制拉格朗日插值函数曲线
Python
注释:支持绘制多条
import matplotlib.pyplot as plt
def function2(x,ltx,j):
cnt=1
for i in range(len(ltx)):
if i==j:
continue
cnt*=(x-ltx[i])
return cnt
def function1(x,ltx,lty):
y=0
for i in range(len(ltx)):
y+=lty[i]/function2(ltx[i],ltx,i)*function2(x,ltx,i)
return y
print("输入注意横纵坐标一一对应。")
ltx=[];lty=[];cnt=0
while True:
cnt+=1;print(f"第{cnt}组即第{cnt}条曲线的输入:")
ltxi=[eval(i) for i in input("输入所有点横坐标,使用空格作为分割:").strip().split()];ltx.append(ltxi)
ltyi=[eval(i) for i in input("输入所有纵横坐标,使用空格作为分割:").strip().split()];lty.append(ltyi)
z=input("继续/结束(Y/N)")
if z=="N":
break
plt.figure(figsize=(16,9))
for i in range(cnt):
x,y=[],[];cnt0=0
while cnt0<=10:
x.append(cnt0)
y.append(function1(cnt0,ltx[i],lty[i]))
cnt0+=0.01
plt.plot(x,y,linewidth=1,label="Curve"+str(i+1))
plt.scatter(ltx[i],lty[i],c="black",s=10)
plt.legend()
plt.grid()
plt.show()
C
注释:1.节点必须最多100个2.由于C无绘图函数所以这里没有绘图。相当于一开始给出n个节点然后就得到了函数,再次输入x就会返回y。支持继续退出输入。
#include <stdio.h>
#define MAX_POINTS 100
double function2(double x,double ltx[],int n,int j) {
double cnt=1;
for (int i=0;i<n;i++){
if(i==j){continue;}
cnt*=(x-ltx[i]);
}
return cnt;
}
double function1(double x,double ltx[],double lty[],int n){
double y=0;
for(int i=0;i<n;i++){
y+=lty[i]/function2(ltx[i],ltx,n,i)*function2(x,ltx,n,i);
}
return y;
}
int main(){
double ltx[MAX_POINTS],lty[MAX_POINTS],x,y;char z;int n;
printf("告诉我有多少节点\n");
scanf("%d",&n);
printf("输入所有点横坐标,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<x[i]);
}
printf("输入所有点纵坐标,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<y[i]);
}
while(1){
printf("请输入x,我将告诉你y:\n");
scanf("%lf",&x);
printf("%lf\n",function1(x,ltx,lty,n));
printf("继续/结束(Y/N)\n");
scanf("%c",&z);
if(z=='N'){break;}
}
}
Matlab
注释:支持绘制多条
%lagelangri
disp("输入注意横纵坐标一一对应。");
ltx=[];lty=[];cnt=0;
while true
cnt=cnt+1;
disp(sprintf("第%d组即第%d条曲线的输入:",cnt,cnt));
ltxi=input("输入所有点横坐标,使用空格作为分割:",'s');ltxi=str2num(ltxi);ltx=[ltx;ltxi];
ltyi=input("输入所有点纵坐标,使用空格作为分割:",'s');ltyi=str2num(ltyi);lty=[lty;ltyi];
z=input("继续/结束(Y/N)", 's');
if z=="N"
break;
end
end
figure('Position',[0,0,1200,900]);
for i=1:cnt
x=[];y=[];cnt0=0;
while cnt0<=10
x=[x,cnt0];y=[y,function1(cnt0,ltx(i,:),lty(i,:))];
cnt0=cnt0+0.01;
end
plot(x,y,'LineWidth',1,'DisplayName',"Curve"+num2str(i));hold on;
scatter(ltx(i,:),lty(i,:),'MarkerFaceColor','black','MarkerEdgeColor','black');hold on;
end
legend('show');grid on;hold off;
%function1
function y=function1(x,ltx,lty)
y=0;
for i=1:length(ltx)
y=y+lty(i)/function2(ltx(i),ltx,i)*function2(x,ltx,i);
end
end
%function2
function cnt=function2(x,ltx,j)
cnt=1;
for i=1:length(ltx)
if i==j
continue
end
cnt=cnt*(x-ltx(i));
end
end
结果
输入
输出
2.3均差与牛顿插值多项式
2.3.1插值多项式的逐次生成
2.3.2均差及其性质
定义2:k阶均差 均差的三性质
2.3.3牛顿插值多项式
牛顿均差插值多项式
例4
这是求均差的例题
注释:我这个比书上的好。因为书上的每一步都有舍入误差。这个是通用的。
Python
ltx=[eval(i) for i in input("输入n个节点的x,使用空格分割").strip().split()]
lt=[[eval(i) for i in input("输入n个节点的y,使用空格分割").strip().split()]]
for i in range(len(ltx)-1):
lt_new=[]
for j in range(len(ltx)-(i+1)):
lt_new.append((lt[i][j+1]-lt[i][j])/(ltx[i+1+j]-ltx[j]))
lt.append(lt_new)
for i in range(len(lt)):
print(f"{i:02d}阶均差:",end=" ")
for j in lt[i]:
print("{:.5f}".format(j),end=" ")
print("")
C
注释:节点必须最多100个
#include <stdio.h>
#define MAX_N 100
double ltx[MAX_N],lt[MAX_N][MAX_N];
int main(){
int n;
printf("输入节点的个数n:");scanf("%d",&n);
printf("输入n个节点的x,使用空格分割:");
for (int i=0;i<n;i++){
scanf("%lf",<x[i]);
}
printf("输入n个节点的y,使用空格分割:");
for (int i=0;i<n;i++){
scanf("%lf",<[0][i]);
}
for (int i=0;i<n-1;i++){
for (int j=0;j<n-i+1;j++){
lt[i+1][j]=(lt[i][j+1]-lt[i][j])/(ltx[i+1+j]-ltx[j]);
}
}
for (int i=0;i<n;i++){
printf("%02d阶均差: ",i);
for (int j=0;j<n-i;j++){
printf("%.5f ",lt[i][j]);
}
printf("\n");
}
return 0;
}
Matlab
ltx=input("输入n个节点的x,使用空格分割:",'s');
ltx=str2num(ltx);
lty=input("输入n个节点的y,使用空格分割:",'s');
lty=str2num(lty);
lt{1}=lty;
for i=1:(length(ltx)-1)
lt_new=zeros(1,length(ltx)-i);
for j=1:(length(ltx)-i)
lt_new(j)=(lt{i}(j+1)-lt{i}(j))/(ltx(i+j)-ltx(j));
end
lt{i+1}=lt_new;
end
for i=1:length(lt)
fprintf("%02d阶均差:",i-1);
fprintf("%.5f ",lt{i});
fprintf("\n");
end
结果
2.3.4差分形式的牛顿插值公式
补充3:牛顿插值
题面描述:任意给出n个节点,求出牛顿插值函数并且绘制牛顿插值函数曲线.
Python
注释:不能支持绘制多条函数曲线
import matplotlib.pyplot as plt
def func2(x,i):
cnt=1
for j in range(i):
cnt*=(x-ltx[j])
return cnt
def func1(x):
cnt=0
for i in range(len(lt_new)):
cnt+=lt_new[i]*func2(x,i)
return cnt
ltx=[eval(i) for i in input("输入n个节点的x,使用空格分割").strip().split()]
lt=[[eval(i) for i in input("输入n个节点的y,使用空格分割").strip().split()]]
for i in range(len(ltx)-1):
lt_new=[]
for j in range(len(ltx)-(i+1)):
lt_new.append((lt[i][j+1]-lt[i][j])/(ltx[i+1+j]-ltx[j]))
lt.append(lt_new)
lt_new=[i[0] for i in lt]
x=[0+i*0.001 for i in range(10001)]
y=[func1(i) for i in x]
plt.figure(figsize=(16,9))
plt.plot(x,y,linewidth=5,color="red")
plt.scatter(ltx,lt[0],s=100,c='blue')
plt.grid();plt.show()
C
注释:1.节点必须最多100个2.由于C无绘图函数所以这里没有绘图。相当于一开始给出n个节点然后就得到了函数,再次输入x就会返回y。支持继续退出输入。
#include <stdio.h>
#define N 100
double ltx[N],lt[N][N];
double func2(double x,int i){
double cnt=1;
for (int j=0;j<i;j++){
cnt*=(x-ltx[j]);
}
return cnt;
}
double func1(double x,int n){
double cnt=0;
for (int i=0;i<n;i++){
double tmp=lt[i][0]*func2(x,i);
cnt+=tmp;
}
return cnt;
}
int main(){
double x;int n;char c;
printf("告诉我有多少节点\n");
scanf("%d",&n);
printf("输入所有点横坐标,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<x[i]);
}
printf("输入所有点纵坐标,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<[0][i]);
}
for(int i=0;i<n-1;i++){
for(int j=0;j<n-1-i;j++){
lt[i+1][j]=(lt[i][j+1]-lt[i][j])/(ltx[i+j+1]-ltx[j]);
}
}
while(1){
printf("输入x,输出y:\n");scanf("%lf",&x);
printf("%lf\n",func1(x,n));
printf("继续/退出(Y/N):\n");char c;scanf(" %c", &c);
if (c=='N'){break;}
}
return 0;
}
Matlab
注释:不能支持绘制多条函数曲线
%niudunchazhi
ltx=input("输入n个节点的x,使用空格分割:",'s');
ltx=str2num(ltx);
lt=input("输入n个节点的y,使用空格分割:",'s');
lt=str2num(lt);
for i=1:length(ltx)-1
lt_new=zeros(1,length(ltx));
for j=1:length(ltx)-i
lt_new(j)=(lt(i,j+1)-lt(i,j))/(ltx(j+i)-ltx(j));
end
lt=[lt;lt_new];
end
lt_new=lt(:,1);
x=0:0.001:10;
y=arrayfun(@(i)func1(i,ltx,lt_new),x);
figure();
plot(x,y,'LineWidth',5,'Color','red');
hold on;
scatter(ltx,lt(1,:),100,'blue','filled');
grid on;
hold off;
%func1
function [cnt]=func1(x,ltx,lt_new)
cnt=0;
for i=1:length(lt_new)
cnt=cnt+lt_new(i)*func2(x,i,ltx);
end
end
%func2
function [cnt]=func2(x,i,ltx)
cnt=1;
for j=1:i-1
cnt=cnt*(x-ltx(j));
end
end
结果
输入
输出
2.4埃尔米特插值
埋坑后面再来解决
2.5分段低次插值
2.5.1高次插值的病态性质
高次插值就是有病,主观上也很好理解
复刻:观察高次插值病态
只有Python因为本质是做可视化并不涉及新的算法(这里算法前面已经提了))
Python
import matplotlib.pyplot as plt
def function12(x,ltx,j):
cnt=1
for i in range(len(ltx)):
if i==j:
continue
cnt*=(x-ltx[i])
return cnt
def function11(x,ltx,lty):
y=0
for i in range(len(ltx)):
y+=lty[i]/function12(ltx[i],ltx,i)*function12(x,ltx,i)
return y
def function02(i):
return 1/(1+pow(i,2))
def function01(n):
lt_new=[-5+10*i/n for i in range(n+1)]
return lt_new
lt=[]
for i in range(2,22,2):
lt.append(function01(i))
for i in lt:
ltx=i;lty=[function02(i) for i in ltx];x=5-5/len(i)
T=function02(x);F=function11(x,ltx,lty);R=T-F
print(f"n={len(i):02d}: {T:.6f} {F:.6f} {R:.6f}")
plt.figure(figsize=(16,9))
ltx=lt[4];lty=[function02(i) for i in ltx]
x,y1,y2=[],[],[];cnt0=-5
while cnt0<=5:
x.append(cnt0)
y1.append(function11(cnt0,ltx,lty))
y2.append(function02(cnt0))
cnt0+=0.01
plt.plot(x,y1,linewidth=1,label="F")
plt.plot(x,y2,linewidth=1,label="T")
plt.scatter(ltx,lty,c="black",s=10)
plt.legend();plt.grid();plt.show()
结果
2.5.2分段线性插值
分段线性插值函数
补充4:分段线性插值
题面描述:任意给出n个节点,求出分段线性插值函数并且绘制分段线性插值函数曲线
Python
注释:支持绘制多条
import matplotlib.pyplot as plt
def function1(x,ltx,lty):
for i in range(len(ltx)):
if ltx[i]>x:
break
return (x-ltx[i])/(ltx[i-1]-ltx[i])*lty[i-1]+(x-ltx[i-1])/(ltx[i]-ltx[i-1])*lty[i]
def function02(i):
return 1/(1+pow(i,2))
def function01(n):
lt_new=[-5+10*i/n for i in range(n+1)]
return lt_new
lt=[]
for i in range(5,56,10):
lt.append(function01(i))
plt.figure(figsize=(16,9))
for i in lt:
ltx=i;lty=[function02(i) for i in ltx]
x,y=[],[];cnt0=-5
while cnt0<=5:
x.append(cnt0)
y.append(function1(cnt0,ltx,lty))
cnt0+=0.01
plt.plot(x,y,linewidth=1,label=str(len(ltx)-1))
ltx=[-5+i*0.01 for i in range(1001)]
lty=[function02(i) for i in ltx]
plt.plot(ltx,lty,linewidth=1,label="TRUE")
plt.legend();plt.grid();plt.show()
C
注释:1.节点必须最多100个2.由于C无绘图函数所以这里没有绘图。相当于一开始给出n个节点然后就得到了函数,再次输入x就会返回y。支持继续退出输入。
#include <stdio.h>
float function1(double x,double ltx[],double lty[],int n){
int i=0;
for (;i<n;i++){
if (ltx[i]>x){
break;
}
}
return (x-ltx[i])/(ltx[i-1]-ltx[i])*lty[i-1]+(x-ltx[i-1])/(ltx[i]-ltx[i-1])*lty[i];
}
int main(){
int n;double x,ltx[100],lty[100];
printf("告诉我有多少节点:\n");scanf("%d",&n);
printf("输入n个节点的x,使用空格分割:\n");
for (int i=0;i<n;i++){
scanf("%lf",<x[i]);
}
printf("输入n个节点的y,使用空格分割:\n");
for (int i=0;i<n;i++) {
scanf("%lf",<y[i]);
}
while(1){
printf("输入x,输出y:\n");scanf("%lf",&x);
printf("%lf\n",function1(x,ltx,lty,n));
printf("继续/退出(Y/N):\n");char c;scanf(" %c", &c);
if (c=='N'){
break;
}
}
return 0;
}
Matlab
%fenduanxianxingchazhi
lt={};
for i=5:10:55
lt{(i-5)/10+1}=function01(i);
end
figure('Position',[0 0 1600 900]);hold on;
for i=1:length(lt)
ltx=lt{i};lty=arrayfun(@function02,ltx);
x=[];y=[];cnt0=-5;
while cnt0<=5
x=[x,cnt0];
y=[y,function3(cnt0,ltx,lty)];
cnt0=cnt0+0.01;
end
plot(x,y,'LineWidth',1,'DisplayName',num2str(length(ltx)-1));
end
ltx=linspace(-5,5,1001);
lty=arrayfun(@function02, ltx);
plot(ltx,lty,'LineWidth',1,'DisplayName','TRUE');
legend;grid on;hold off;
%function01
function lt_new=function01(n)
lt_new=linspace(-5,5,n+1);
end
%function02
function result=function02(i)
result=1/(1+i^2);
end
%function3
function result=function3(x,ltx,lty)
for i=1:length(ltx)
if ltx(i) > x
break;
end
end
result=(x-ltx(i))/(ltx(i-1)-ltx(i))*lty(i-1)+(x-ltx(i-1))/(ltx(i)-ltx(i-1))*lty(i);
end
结果
这个不够光滑引出分段三次埃尔米特插值
2.5.3分段三次埃尔米特插值
补充5:分段三次埃尔米特插值
题面描述:任意给出n个节点,求出分段三次埃尔米特插值函数并且绘制分段三次埃尔米特插值函数曲线。注意输入的x与y与导数值一一对应。
Python
注释:不能支持绘制多条函数曲线
import matplotlib.pyplot as plt
def function_this(x):
for i in range(len(ltx)):
if ltx[i]>x:
break
cnt1=(x-ltx[i])/(ltx[i-1]-ltx[i])
cnt2=(x-ltx[i-1])/(ltx[i]-ltx[i-1])
return pow(cnt1,2)*(1+2*cnt2)*lty[i-1]+pow(cnt2,2)*(1+2*cnt1)*lty[i]+\
pow(cnt1,2)*(x-ltx[i-1])*ltz[i-1]+pow(cnt2,2)*(x-ltx[i])*ltz[i]
ltx=[eval(i) for i in input("输入所有点横坐标,使用空格作为分割:").strip().split()]
lty=[eval(i) for i in input("输入所有纵横坐标,使用空格作为分割:").strip().split()]
ltz=[eval(i) for i in input("输入所有点导数值,使用空格作为分割:").strip().split()]
plt.figure(figsize=(16,9))
x=[-5+i*0.01 for i in range(1001)]
y=[function_this(i) for i in x]
plt.plot(x,y,linewidth=2,color="b")
plt.scatter(ltx,lty,s=20,c="r")
plt.grid();plt.show()
分段三次埃尔米特插值C
注释:1.节点必须最多100个2.由于C无绘图函数所以这里没有绘图。相当于一开始给出n个节点然后就得到了函数,再次输入x就会返回y。支持继续退出输入。
#include<stdio.h>
#include<math.h>
double ltx[100],lty[100],ltz[100];int n;
double function_this(double x){
int i;double cnt1,cnt2;
for(i=0;i<n;i++){
if(ltx[i]>x){
break;
}
}
cnt1=(x-ltx[i])/(ltx[i-1]-ltx[i]);
cnt2=(x-ltx[i-1])/(ltx[i]-ltx[i-1]);
return pow(cnt1,2)*(1+2*cnt2)*lty[i-1]+pow(cnt2,2)*(1+2*cnt1)*lty[i]+
pow(cnt1,2)*(x-ltx[i-1])*ltz[i-1]+pow(cnt2,2)*(x-ltx[i])*ltz[i];
}
int main()
{
printf("告诉节点个数:\n");scanf("%d",&n);
printf("输入所有点横坐标,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<x[i]);
}
printf("输入所有点纵坐标,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<y[i]);
}
printf("输入所有点导数值,使用空格作为分割:\n");
for(int i=0;i<n;i++){
scanf("%lf",<z[i]);
}
while(1)
{
double x;char z;
printf("输入x输出y\n");scanf("%lf",&x);
printf("%lf\n",function_this(x));
printf("继续/结束(Y/N)\n");scanf("%c",&z);
if(z=='N'){
break;
}
}
return 0;
}
Matlab
注释:不能支持绘制多条函数曲线
%fenduansanciaiermitechazhi
ltx=input("输入所有点横坐标,使用空格作为分割:",'s');
ltx=str2num(ltx);
lty=input("输入所有点纵坐标,使用空格作为分割:",'s');
lty=str2num(lty);
ltz=input("输入所有点导数值,使用空格作为分割:",'s');
ltz=str2num(ltz);
x=-5:0.01:5;
y=arrayfun(@(x)function_this(x,ltx,lty,ltz),x);
figure('Position',[0 0 800 500]);
plot(x,y,'LineWidth',2,'Color','b');
hold on;
scatter(ltx,lty,20,'r','filled');
grid on;
%function_this
function output=function_this(x,ltx,lty,ltz)
for i=1:length(ltx)
if ltx(i)>x
break
end
end
cnt1=(x-ltx(i))/(ltx(i-1)-ltx(i));
cnt2=(x-ltx(i-1))/(ltx(i)-ltx(i-1));
output=cnt1^2*(1+2*cnt2)*lty(i-1)+cnt2^2*(1+2*cnt1)*lty(i)+...
cnt1^2*(x-ltx(i-1))*ltz(i-1)+cnt2^2*(x-ltx(i))*ltz(i);
end
结果
输入
输出
同样不够光滑引出三次样条插值
2.6三次样条插值
直接跳过,根本不会
评注
复习与思考题
习题
计算实习题
计算实习题1
代码也就不给了吧,其实就是前面牛顿插值法改改就好了,结果如下(跳过三次样条插值):
计算实习题2
代码也就不给了吧,其实就是前面多项式插值法改改就好,结果如下(跳过三次样条插值):
计算实习题3
代码也就不给了吧,其实就是前面拉格朗日插值法改改就好了,结果如下(跳过三次样条插值):