插值算法Word文件下载.docx
《插值算法Word文件下载.docx》由会员分享，可在线阅读，更多相关《插值算法Word文件下载.docx（26页珍藏版）》请在冰豆网上搜索。
不少实际插值问题不仅要求函数值相等，而且还要求导数值也相等。
这就导致下面的Hermite插值。
给定
（4）三角函数插值
当被插函数是以2π为周期的函数时，通常用n阶三角多项式作为插值函数，并通过高斯三角插值表出。
（5）其它提法
插值（Interpolation），有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。
有些相机使用插值，人为地增加图像的分辨率。
插值：
用来填充图像变换时像素之间的空隙。
（6）三次样条插值:
上面介绍的分段线性插值,其总体光滑程度不够.在数学上,光滑程度的定量描述是：
函数（曲线）的k阶导数存在且连续,则称该曲线具有k阶光滑性.自然,阶数越高光滑程度越好.于是,分段线性插值具有零阶光滑性,也就是不光滑；
分段三次埃尔米特插值具有一阶光滑性.仅有这些光滑程度,在工程设计和机械加工等实际中是不够的.提高分段函数如多项式函数的次数,可望提高整体曲线的光滑程度.但是,是否存在较低次多项式达到较高阶光滑性的方法？
三次样条插值就是一个很好的例子.
样条曲线本身就来源于飞机、船舶等外形曲线设计中所用的绘图工具.在工程实际中,要 求这样的曲线应该具有连续的曲率,也就是连续的二阶导数.值得注意的是分段插值曲线的光滑性关键在于段与段之间的衔接点（节点）处的光滑性.
三次样条函数记为s（x）,它是定义在区间[a,b]上的函数,满足:
1）s（x）在每一个小区间[Xi-1,Xi]上是一个三次多项式函数；
2）在整个区间[a,b]上,其二阶导数存在且连续.即在每个节点处的二阶导数连续.
三次样条插值问题的提法:
给定函数f（x）在n+1个节点x0,x1,...,xn处的函数值为y0,y1,...,yn,求一个三次样条函数s（x）,使其满足：
s（xi）=yi,i=0,1,…,n.
如何确定三次样条函数在每一个小区间上的三次多项式函数的系数呢？
这是一个比较复杂的问题,这里只介绍确定系数的思想.
分段线性插值在每一段的线性函数的两个参数,是由两个方程（两个端点处的函数值为给定值）唯一确定；
对于三次样条插值呢,每一个区间上的三次函数有四个参数,而在该区间上由两个端点的函数值为给定值只能够产生两个方程,仅此不足以唯一确定四个参数.注意到三次样条函数对整体光滑性的要求,其二阶导数存在且连续,从全局的角度上考虑参数个数与方 程个数的关系如下：
参数：
每个小段上4个,n个小段共计4n个.
方程：
每个小段上由给定函数值得到2个,n个小段共计2n个；
光滑性要求每一个内部节点处的一阶、二阶导数连续,得出其左右导数相等,因此,每个节点产生2个方程,共计2（n-1）个.
现在得到了4n-2个方程,还差两个.为此,常用的方法是对边界节点除函数值外附加要求.这就是所谓的边界条件.需要两个,正好左右两个端点各一个.常用如下三类边界条件.m边界条件：
s'
（X0）=m0,s'
（Xn）=mn.即两个边界节点的一阶导数值为给定值：
m0,mn.
M边界条件：
'
（x0）=m0,s'
（xn）=mn.即两个边界节点的二阶导数值为给定值：
m0,mn.特别地,当m0和mn都为零时,称为自然边界条件.
周期性边界条件：
（x0）=s'
（xn）;
s'
（xn）.
以上分析说明,理论上三次样条插值函数是确定的,具体如何操作,可以查阅有关资料.
例题，观测数据为
x=[012345678910];
y=[020-4040-2031];
待求的三次多项式函数s（x）在[010]上有连续的一阶，二阶导数.我们通过简单的讨论来认识问题。
在第一区间[01]、第二区间[12]上考虑两个三次多项式
s（x）=s1*x^3+s2*x^2+s3*x+s4
r（x）=r1*x^3+r2*x^2+r3*x+r4
示意图:
可以得到
s（0）=s1*0^3+s2*0^2+s3*0+s4=0
（1）
s
（1）=s1*1^3+s2*1^2+s3*1+s4=2
（2）
r
（1）=r1*2^3+r2*2^2+r3*2+r4=2（3）
r
（2）=r1*2^3+r2*2^2+r3*2+r4=0（4）
一阶导函数
（x）=3*s1*x^2+2*s2*x+s3
r'
（x）=3*r1*x^2+2*r2*x+r3
由一阶导数的连续性且在1点处相等，有
3*s1*1^2+2*s2*1+s3=3*r1*1^2+2*r2*1+r3（5）
二阶导函数
（x）=6*s1*x+2*s2
（x）=6*r1*x+2*r2
由二阶导数的连续性且在1点处相等，有
6*s1*1+2*s2=6*r1*1+2*r2（6）
由m边界条件
（0）=1.6,r'
（2）=0.3
有
3*s1*0^2+2*s2*0+s3=1.6（7）
3*r1*2^2+2*r2*2+r3=0.3（8）
M边界条件
（0）=-1,r'
（2）=1
有
6*s1*0+2*s2=-1（7'
）
6*r1*2+2*r2=1（8'
由周期性边界条件
（0）=r'
（2）
3*s1*0^2+2*s2*0+s3=-1（7'
）
3*r1*2^2+2*r2*2+r3=1（8'
）这样，对于两个多项式的8个未知量，我们给出了8个方程。
三次样条曲线的难点在于，我们不能分段去求解方程，完成绘图。
//三次样条插值函数代码：
/*
注:
所用的数据表如下:
x0701302103375787761012114214621841
y05778103135182214244256272275
边界条件为:
（1）y'
（0）=1,y'
（1841）=0;
（2）y"
（0）=0,y"
（1841）=0.
求解在以上边界条件下在x[k]=50k（k=1,2...,36）处的函数值程序源代码*/
二:
插值算法的代码实现：
1三次样条的代码实现及测试结果
//求三次样条插值函数
#include<
iostream>
iomanip>
usingnamespacestd;
constintMAX=50;
floatx[MAX],y[MAX],h[MAX];
//变量设置：
x为各点横坐标；
y为各点纵坐标；
h为步长
floatc[MAX],a[MAX],fxym[MAX];
floatf（intx1,intx2,intx3）/*****************求差分函数（含三个参数）****************************/
{
floata=（y[x3]-y[x2]）/（x[x3]-x[x2]）;
floatb=（y[x2]-y[x1]）/（x[x2]-x[x1]）;
return（a-b）/（x[x3]-x[x1]）;
}
voidcal_m（intn）/***********************用追赶法求解出弯矩向量M……***************************/
{
floatB[MAX];
B[0]=c[0]/2;
for（inti=1;
i<
n;
i++）
B[i]=c[i]/（2-a[i]*B[i-1]）;
//fxym[0]=fxym[0]/2;
=n;
fxym[i]=（fxym[i]-a[i]*fxym[i-1]）/（2-a[i]*B[i-1]）;
for（inti=n-1;
i>
=0;
i--）
fxym[i]=fxym[i]-B[i]*fxym[i+1];
}
voidprintout（intn）;
intmain（）
intn,i;
charch;
do
{
cout<
<
"
请输入已知断点个数:
;
cin>
>
n;
for（i=0;
{
cout<
PleaseputinX"
i<
:
cin>
x[i];
//cout<
endl;
PleaseputinY"
y[i];
}
for（i=0;
i++）//求步长；
其数组值较之x，y个数少一
h[i]=x[i+1]-x[i];
Please输入边界条件\n1:
已知两端的一阶导数\n2:
两端的二阶导数已知\n默认:
自然边界条件\n"
intt;
floatf0,f1;
t;
switch（t）
case1:
cout<
PleaseputinY0\'
Y"
n<
\'
\n"
//显示数据为Y0'
至Yn'
，即断点的一阶导数
f0>
f1;
c[0]=1;
a[n]=1;
fxym[0]=6*（（y[1]-y[0]）/（x[1]-x[0]）-f0）/h[0];
fxym[n]=6*（f1-（y[n]-y[n-1]）/（x[n]-x[n-1]））/h[n-1];
break;
case2:
PleaseputinY0\"
\"
//显示数据为Y0"
至Yn"
，即断点的二阶导数
c[0]=a[n]=0;
fxym[0]=2*f0;
fxym[n]=2*f1;
default:
不可用\n"
//待定
};
//switch
for（i=1;
fxym[i]=6*f（i-1,i,i+1）;
//调用差分函数（only!
a[i]=h[i-1]/（h[i]+h[i-1]）;
c[i]=1-a[i];
a[n]=h[n-1]/（h[n-1]+h[n]）;
cal_m（n）;
//调用弯矩函数（only！
\n输出三次样条插值函数：
printout（n）;
//调用求解三次样条插值函数；
函数输出
Doyoutohaveanthertry?
y/n:
ch;
}
while（ch=='
y'
||ch=='
Y'
）;
return0;
voidprintout（intn）/***************求三次样条插值函数（因已知断点个数而异）***********************/
cout<
setprecision（6）;
//通过操作器setprecision（）设置有效位数；
其为头文件<
iomanip.h>
所包含；
括号内为参数。
for（inti=0;
i++）//所输出函数个数由所设断点个数而定
{
i+1<
["
x[i]<
"
x[i+1]<
]\n"
\t"
floatt=fxym[i]/（6*h[i]）;
if（t>
0）
t<
*（"
-x）^3"
else
-t<
*（x-"
）^3"
t=fxym[i+1]/（6*h[i]）;
+"
-"
\n\t"
t=（y[i]-fxym[i]*h[i]*h[i]/6）/h[i];
+"
-x）"
else
-"
t=（y[i+1]-fxym[i+1]*h[i]*h[i]/6）/h[i];
）"
endl<
}
/****************提供测试数据：
（来自课本页例《数值分析》清华大学出版社第四版）***********/
/*输入
3
27.74.1
284.3
294.1
303.0
1
3.0-4.0*/
/*输出
输出三次样条插值函数：
[27.7,28]
13.07*（x-28）^3+0.22*（x-27.7）^3
+14.84*（28-x）+14.31*（x-27.7）
[28,29]
0.066*（29-x）^3+0.1383*（x-28）^3
+4.234*（29-x）+3.962*（x-28）
3:
[29,30]
0.1383*（30-x）^3-1.519*（x-29）^3
+3.962*（30-x）+4.519*（x-29）
*/
/*利用三次样条插值法求曲线在某个点处的函数值（本题用第二个边界条件）*/
＃include<
stdio.h>
math.h>
#defineN11//N表示节点的个数
voidmain（）
doublex[N]={0,70,130,210,337,578,776,1012,1142,1462,1841},
y[N]={0,57,78,103,135,182,214,244,256,272,275},
h[N],a[N],b[N],A[N],B[N],m[N],s[37],xx[37];
//s[N]表示曲线,其中数组xx[N]表示自变量x[k]
inti,k;
h[0]=x[1]-x[0];
//初始化h0,a0,b0,A0,B0
a[0]=1;
b[0]=3*（y[1]-y[0]）/h[0];
A[0]=-a[0]/2;
B[0]=b[0]/2;
for（i=0;
N;
i++）//求hi
h[i]=x[i+1]-x[i];
for（i=1;
N-1;
i++）{//求ai,bi
a[i]=h[i-1]/（h[i-1]+h[i]）;
b[i]=3*（（1-a[i]）*（y[i]-y[i-1]）/h[i-1]+a[i]*（y[i+1]-y[i]）/h[i]）;
i++）{//求Ai,Bi
A[i]=-a[i]/（2+（1-a[i]）*A[i-1]）;
B[i]=（b[i]-（1-a[i]）*B[i-1]）/（2+（1-a[i]）*A[i-1]）;
m[N-1]=（b[N-1]-（1-a[N-1]）*B[N-2]）/（2+（1-a[N-1]）*A[N-2]）;
//求Mn的值
for（i=N-2;
i>
=0;
i--）//求m0,m1,-----mn-1的值
m[i]=A[i]*m[i+1]+B[i];
for（k=1;
k<
=36;
k++）{//求曲线在x[k]处的函数值
xx[k]=50*k;
for（i=0;
i++）//找出x[k]所在的区间
if（xx[k]>
=x[i]&
&
xx[k]<
=x[i+1]）
//s[k]即为x[k]所在区间的三次样条插值函数,以下即为求在x[k]处的函数值
s[k]=（1+2*（xx[k]-x[i]）/（x[i+1]-x[i]））*pow（（xx[k]-x[i+1]）/（x[i]-x[i+1]）,2）*y[i]+
（1+2*（xx[k]-x[i+1]）/（x[i]-x[i+1]））*pow（（xx[k]-x[i]）/（x[i+1]-x[i]）,2）*y[i+1]+
（xx[k]-x[i]）*pow（（xx[k]-x[i+1]）/（x[i]-x[i+1]）,2）*m[i]+
（xx[k]-x[i+1]）*pow（（xx[k]-x[i]）/（x[i+1]-x[i]）,2）*m[i+1];
printf（"
所求的外形曲线在x[k]=50*k（k=1,2,...,36）处的函数值分别为:
k++）//输出在x[k]处的函数值
printf（"
s（%4d）=%13.8f\n"
50*k,s[k]）;
printf（"
m[i]的值分别为:
i++）
m（%2d）=%8f\n"
i,m[i]）;
//程序代码实现2：
#definen4//三次样条插值分段数=数据的组数减一
voidchase（doublea[n-1][n-1],doublef[n-1],doublem[n+1]）
//三对角阵的Crout分解
doubleL[n-1][n-1];
doubleU[n-1][n-1];
doubley[n-1];
doublex[n-1];
intt,r=0;
for（inti=0;
n-1;
for（intj=0;
j<
j++）
L[i][j]=U[i][j]=0;
i++）//判断是否符合Crout分解定理
if（i==0）
if（a[i][i]==0||a[i][i+1]==0）
{
cout<
a["
n-1<
]["
]可约"
break;
}
if（fabs（a[i][i]）<
=fabs（a[i][i+1]））
]强超"
break;
elseif（i==n-1-1）
if（a[i][i]==0||a[i][i-1]==0）
=fabs（a[i][i-1]））
if（a[i][i]==0||a[i][i-1]==0||a[i][i+1]==0）
cout
升级会员 