数据插值和函数逼近
1 数据插值
由已知样本点，以数据更为平滑为目标，求出其他点处的函数
值。在信号处理与图像处理上应用广泛。
求解方法：
y1=interp1(x,y,x1,'方法')
z1=interp2(x,y,z,x1,y1,'方法')
1.1 一维数据的插值
例：假设样本点来自x e x x x f x sin )53()(52-+-=，进行插值处理，得到平
例：草图样条曲线功能。
function sketch()
x=[]; y=[]; gca;
hold on; axis([0 1,0 1]);
while 1
[x0,y0,button]=ginput(1);
if(isempty(button)) break; end;
x=[x x0]; y=[y y0];
plot(x,y,'*');
end;
xx=[x(1):(x(end)-x(1))/100:x(end)];
yy=interp1(x,y,xx,'spline');
plot(xx,yy,x,y,'*');
end
1.2 二维网格数据的插值
例3：假设样本点来自xy y x e x x z ----=2
2)2(2，进行插值处理，得到平滑的曲
1.3 二维一般分布数据的插值
例4：假设样本点来自xy y x e
x x z ----=22)2(2，进行插值处理，得到平滑的曲
2 样条插值函数逼近
由已知样本点，求能对其较好拟合的函数表达式。
求解方法：
S=csapi(x,y); % 定义一个三次样条函数类
S=spapi(k,x,y); % 定义一个k 次B 样条函数类
ys=fnval(S,xs); % 计算插值结果
fnplt(S); % 绘制插值结果
例：从)sin(x y =中取样本点，计算三次样条函数。
练习：利用上例中的变量S，生成各个曲线段的参数方程，并绘
例：从)sin(x y =中取样本点，计算三次B 样条函数。
3 多项式的插值与函数逼近
3.1 多项式的基本运算
设多项式为：116251234++-x x x 直接输入系数向量创建多项式
p1=[1 -12 0 25 116];
y=poly2sym(p1);
求多项式的根
r=roots(p1);
由多项式的根逆推多项式
p2=poly(r); p2和p1等价
求多项式的值
polyval (p1,[1 2; 3 4]); % 以数组为运算单位
polyvalm(p1,[1 2; 3 4]); % 以矩阵为运算单位
求多项式的微积分
q=polyder(p); % 导数
r=polyint(q); % 积分
设多项式为：322++x x 与2+x 。
a=[1 2 3]; b=[1 2];
c=conv(a,b); % 乘积
[q,r]=deconv(a,b); % 相除
q=1 0 ％商多项式
r=0 0 3 ％余多项式
3.2 Lagrange 插值算法
∑∏=≠=--=N i N i j j
i i
x x x x y x 1,1)(φ
3.3 多项式函数逼近
求解方法：
p=polyfit (x,y,n); 计算n 次多项式的降幂排列的系数向量 y=polyval (p,x); 计算多项式的值
主要用于针对实验数据，进行解析描述。
例：假设样本点来自x e x x x f x sin )53()(52-+-=，比较3次、5次多项式插
4 函数设计中的参数逼近
4.1 函数线性组合的参数逼近
已知某函数是若干已知函数的组合，
C F x f c x f c x f c y n n ?=+++=)()()(2211 ，已知数据
),(11y x ,),(22y x ,……),(M M y x ，则有方程组：
Y C F =?
????????????=?????
???????=????????????=M n M n M M n n y y y Y c c c C x f x f x f x f x f x f x f x f x f F 2121212222
111211,,)()()()()()()()()( 有n 个待求的未知系数，有M 个方程，所以只能求方程的最小
二乘解： Y F C \=
例：已知函数原型为2443321)2cos()(x c e x c e
c c x g x x +-++=--，已知以下数
C
4.2 函数中的参数逼近
已知数据),(11y x ,),(22y x ,……),(N N y x ，满足函数),(x a f ，其中a
是待定系数的向量。目标是求出a 。
即定义最优化问题：∑=-=N
i i i a x a f y J 12)],([min
[a,res]=lsqcurvefit(fun,a0,x,y)
例：设采样点如下：
x=0:1:10;
y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
已知数据满足)sin(*),(5314
2x a e a e a x a f x a x a --+=，求最佳拟合的参
数。
4 多项式运算
数据插值和函数逼近 1 数据插值 由已知样本点，以数据更为平滑为目标，求出其他点处的函数 值。在信号处理与图像处理上应用广泛。 求解方法： y1=interp1(x,y,x1,'方法') z1=interp2(x,y,z,x1,y1,'方法') 1.1 一维数据的插值 例：假设样本点来自x e x x x f x sin )53()(52-+-=，进行插值处理，得到平 例：草图样条曲线功能。 function sketch() x=[]; y=[]; gca; hold on; axis([0 1,0 1]); while 1 [x0,y0,button]=ginput(1);
if(isempty(button)) break; end; x=[x x0]; y=[y y0]; plot(x,y,'*'); end; xx=[x(1):(x(end)-x(1))/100:x(end)]; yy=interp1(x,y,xx,'spline'); plot(xx,yy,x,y,'*'); end 1.2 二维网格数据的插值 例3：假设样本点来自xy y x e x x z ----=2 2)2(2，进行插值处理，得到平滑的曲
1.3 二维一般分布数据的插值 例4：假设样本点来自xy y x e x x z ----=22)2(2，进行插值处理，得到平滑的曲 2 样条插值函数逼近 由已知样本点，求能对其较好拟合的函数表达式。 求解方法： S=csapi(x,y); % 定义一个三次样条函数类 S=spapi(k,x,y); % 定义一个k 次B 样条函数类 ys=fnval(S,xs); % 计算插值结果 fnplt(S); % 绘制插值结果 例：从)sin(x y =中取样本点，计算三次样条函数。
计算方法实验—三次样条插值 机电学院075094-19 苏建加 20091002764 题目：求压紧三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且一阶导 数边界条件S'(-3)=-1;S'（4）=1。 解：首先计算下面的值： 记 1--=j j j x x h ； 1++=j j j j h h h u ；1=+j j u λ ； ?? ????????---+=-++++-j j j j j j j j j j j h y y h y y h h x x x f 1111 111],,[ ；M j =)(''j x s ;],,[611+-=j j j j x x x f d ; h1=-2-(-3)=1;h2=1-(-2)=3;h3=4-1=3; u1=1/4;u2=3/6; d1=6/4*(3/3-(-2)/1)=4.5;d2=6/6*(-2/3-3/3)=-5/3; 由于边界条件S'(-3)=-1;S'（4）=1,得到如下 式子： d0=6/1*(-2/1-(-1))=-6； d3=6/3*(1-(-2)/3)=10/3; 所以得到4个含参数m0~m3 的线性代数方程组为： 2.0000 1.0000 0 0 m0 0.2500 2.0000 0.7500 0 m1 0 0.5000 2.0000 0.5000 m2 0 0 1.0000 2.0000 m3 利用matlab 求解方程得： m = -4.9032 3.8065 -2.5161 2.9247 所以 S1(x)=-0.8172*(-2-x)^3+ 0.6344*(x+3)^3+2.8172*(-2-x)-0.6344*(x+3) x ∈[-3，-2] S2(x)=0.2115*(1-x)^3 -0.1398*(x+2)^3- 1.9032*(1-x)+ 2.2581*(x+2) x ∈[-2，1] S3(x)=-0.1398*(4-x)^3+0.1625(x-1)^3+ 2.2581*(4-x)-1.1290*(x-1) x ∈[1，4] 化简后得：S1(x)=1.4516*x^3 + 10.6128*x^2 + 23.4836*x + 16.1288 x ∈[-3，-2] S2(x)=-0.3513x^3-0.2043x^2+1.8492x+1.7061 x ∈[-2，1] S3(x)=0.3023x^3-2.1651x^2+3.8108x+1.0517 x ∈[1，4] 画图验证：
佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 10 指导教师 陈剑 成 绩 日 期 月 日 一、实验目的 1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法； 2、讨论插值的Runge 现象 3、学会Matlab 提供的插值函数的使用方法，会用这些函数解决实际问题。 二、实验原理 1、拉格朗日插值多项式 2、牛顿插值多项式 3、三次样条插值 三、实验步骤 1、用MATLAB 编写独立的拉格朗日插值多项式函数 2、用MATLAB 编写独立的牛顿插值多项式函数 3、用MATLAB 编写独立的三次样条函数（边界条件为第一、二种情形） 4、已知函数在下列各点的值为： 根据步骤1,2,3编好的程序，试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值，并用图给出 ｛(,),0.20.08,0,1,2, ,10i i i x y x i i =+=｝，4()L x 、4()P x 和()S x 。 5、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数 2 1 (),(11)125f x x x = -≤≤+作多项式插值，对不同n 值，分别画出插值函数及()f x 的图形。 6、下列数据点的插值
可以得到平方根函数的近似，在区间[0,64]上作图。 （1）用这9个点作8次多项式插值8()L x 。 （2）用三次样条（第一边界条件）程序求()S x 。 7、对于给函数2 1 ()125f x x = +在区间[-1,1]上取10.2(0,1, ,10)i x i i =-+=，试求3次 曲线拟合，试画出拟合曲线并打印出方程，与第5题的结果比较。 四、实验过程与结果： 1、Lagrange 插值多项式源代码： function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化 %循环方式求L 系数，并求和： for i = 1:length(y) for j = 1:length(x) if i ~= j mu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end end ya = ya + y(i) * mu ; mu = 1; end 2、Newton 源代码： function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化：
（1）编写三条样条插值函数程序如下： x=[1 4 9 16 25 36 49 64 81]; y=[1 2 3 4 5 6 7 8 9]; n=length(x); lamda(1)=1; miu(n)=1; h=diff(x); df=diff(y)./diff(x); d(1)=6*(df(1)-1/2)/h(1); d(n)=6*(0.5*81^-0.5-df(n-1))/h(n-1); for j=2:n-1 lamda(j)=h(j)/(h(j-1)+h(j)); miu(j)=h(j-1)/(h(j-1)+h(j)); d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j)); end miu=miu(2:end); u=diag(miu,-1);r=diag(lamda,1);a=diag(2*ones(1,n)); A=u+r+a; %求出矩阵形式的线性方程组 M=inv(A)*d'; %求出M值 syms g for j=1:n-1 s(j)=M(j)*(x(j+1)-g)^3/(6*h(j))+M(j+1)*((g-x(j))^3/(6*h(j)))+(y(j)-M( j)*h(j)^2/6)*(x(j+1)-g)/h(j)+(y(j+1)-M(j+1)*h(j)^2/6)*(g-x(j))/h(j); end format rat for j=1:n-1 S(j,:)=sym2poly(s(j)); %三条样条插值函数 end %生成三次样条插值函数图象 for j=1:n-1 x1=x(j):0.01:x(j+1); y1=polyval(S(j,:),x1); plot(x1,y1,x,y,'o'); title('spline 三次样条插值函数图象'); xlabel('x'); ylabel('y'); grid on; hold on; end
MATLAB三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（1，2，…，n）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先，通过函数 sanwanj得到方程的系数矩阵，即追赶法方程的四个向量参数，接下来调用 追赶法（在intersanwj函数中），得到三次样条分段函数系数因子，然后进行多项式合并得 到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算改 点的值。附：追赶法程序 chase %%%%%%%%%%%%%% function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b) % 三弯矩样条插值 % 将插值点分两次输入，x0 y0 单独输入 % 边值条件a的二阶导数 y1a 和b的二阶导数 y1b，这里建议将y1a和y1b换成y2a和 y2b，以便于和三转角代码相区别 n=length(x);m=length(y); if m~=n error('x or y 输入有误，再来'); end v=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1); w=2*ones(n+1); h0=x(1)-x0; h=zeros(n-1,1); for k=1:n-1 h(k)=x(k+1)-x(k); end v(1)=h0/(h0+h(1)); u(1)=1-v(1); d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1)); % for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k)); u(k)=1-v(k); d(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k)); end newv=[v;1]; newu=[1;u]; d0=6*((y(1)-y0)/h0-y1a)/h0;
Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下：
问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得：
MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤
（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x)_____________(1) %第一类边界条件下三次样条插值； %xi所求点； %yi所求点函数值； %x已知插值点； %y已知插值点函数值； %f_0左端点一次导数值； %f_n右端点一次导数值； n = length(x0); z = length(y0); h = zeros(n-1,1); k=zeros(n-2,1); l=zeros(n-2,1); S=2*eye(n); fori=1:n-1 h(i)= x0(i+1)-x0(i); end fori=1:n-2 k(i)= h(i+1)/(h(i+1)+h(i)); l(i)= 1-k(i);
end %对于第一种边界条件： k = [1;k];_______________________(2) l = [l;1];_______________________(3) %构建系数矩阵S： fori = 1:n-1 S(i,i+1) = k(i); S(i+1,i) = l(i); end %建立均差表： F=zeros(n-1,2); fori = 1:n-1 F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i)); end D = zeros(n-2,1); fori = 1:n-2 F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i)); D(i,1) = 6 * F(i,2); end %构建函数D： d0 = 6*(F(1,2)-f_0)/h(1);___________(4)
非常类似前面的三弯矩法，这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj，追赶法程序通用，代码如下。 %%%%%%%%%%%%%%%%%%% function [newu,w,newv,d]=sanzhj(x,y,x0,y0,y1a,y1b) % 三转角样条插值 % 将插值点分两次输入，x0 y0 单独输入 % 边值条件a的一阶导数 y1a 和b的一阶导数 y1b n=length(x);m=length(y); if m~=n error('x or y 输入有误，再来'); end v=ones(n-1,1); u=ones(n-1,1); d=zeros(n-1,1); w=2*ones(n-1,1); h0=x(1)-x0; h=zeros(n-1,1); for k=1:n-1 h(k)=x(k+1)-x(k); end v(1)=h0/(h0+h(1)); u(1)=1-v(1); d(1)=3*(v(1)*(y(2)-y(1))/h(1)+u(1)*((y(1)-y0))/h0); % for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k)); u(k)=1-v(k); d(k)=3*(v(k)*(y(k+1)-y(k))/h(k)+u(k)*(y(k)-y(k-1))/h(k-1)); end d(1)=d(1)-u(1)*y1a; d(n-1)=d(n-1)-v(n-1)*y1b; newv=v(1:n-2,:); newu=u(2:n-1,:); %%%%%%%%%%%% function intersanzhj(x,y,x0,y0,y1a,y1b) % 三转角样条插值
13. 数据插值与拟合 实际中，通常需要处理实验或测量得到的离散数据（点）。插值与拟合方法就是要通过离散数据去确定一个近似函数（曲线或曲面），使其与已知数据有较高的拟合精度。 1.如果要求近似函数经过所已知的所有数据点，此时称为插值问 题（不需要函数表达式）。 2.如果不要求近似函数经过所有数据点，而是要求它能较好地反 映数据变化规律，称为数据拟合（必须有函数表达式）。 插值与拟合都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数。区别是：【插值】不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。【拟合】要求得到一个具体的近似函数的表达式。 因此，当数据量不够，但已知已有数据可信，需要补充数据，此时用【插值】。当数据基本够用，需要寻找因果变量之间的数量关系（推断出表达式），进而对未知的情形作预测，此时用【拟合】。
一、数据插值 根据选用不同类型的插值函数，逼近的效果就不同，一般有：（1）拉格朗日插值（lagrange插值） （2）分段线性插值 （3）Hermite （4）三次样条插值 Matlab 插值函数实现： （1）interp1( ) 一维插值 （2）intep2( ) 二维插值 （3）interp3( ) 三维插值 （4）intern( ) n维插值 1.一维插值（自变量是1维数据） 语法：yi = interp1(x0, y0, xi, ‘method’) 其中，x0, y0为原离散数据（x0为自变量，y0为因变量）；xi为需要插值的节点，method为插值方法。 注：（1）要求x0是单调的，xi不超过x0的范围； （2）插值方法有‘nearest’——最邻近插值；‘linear’——线性插值；‘spline’——三次样条插值；‘cubic’——三次插值；
Matlab中插值函数汇总和使用说明 MATLAB中的插值函数为interp1，其调用格式 为： yi= interp1(x,y,xi,'method') 其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量， 'method'表示采用的插值方法，MATLAB提供的插值方法有几种： 'method'是最邻近插值， 'linear'线性插值； 'spline'三次样条插值； 'cubic'立方插值．缺省时表示线性插值 注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。 例如：在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28，27，25，20，18，15，13， 推测中午12点（即13点）时的温度． x=0:2:24; y=[12 9 9 10 18 24 28 27 25 20 18 15 13]; a=13; y1=interp1(x,y,a,'spline') 结果为： 27.8725 若要得到一天24小时的温度曲线，则： xi=0:1/3600:24; yi=interp1(x,y,xi, 'spline'); plot(x,y,'o' ,xi,yi)
命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。 x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。(2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数pchip，用于对
(第一边界条件源代码: function y=yt1(x0,y0,f_0,f_n,x _____________(1 %第一类边界条件下三次样条插值; %xi 所求点; %yi所求点函数值; %x 已知插值点; %y 已知插值点函数值; %f_0左端点一次导数值; %f_n右端点一次导数值; n = length(x0; z = length(y0; h = zeros(n-1,1; k=zeros(n-2,1; l=zeros(n-2,1; S=2*eye(n; fori=1:n-1 h(i= x0(i+1-x0(i; end fori=1:n-2
k(i= h(i+1/(h(i+1+h(i; l(i= 1-k(i; end %对于第一种边界条件: k = [1;k]; _______________________(2 l = [l;1]; _______________________(3 %构建系数矩阵 S : fori = 1:n-1 S(i,i+1 = k(i; S(i+1,i = l(i; end %建立均差表: F=zeros(n-1,2; fori = 1:n-1 F(i,1 = (y0(i+1-y0(i/(x0(i+1-x0(i; end D = zeros(n-2,1; fori = 1:n-2 F(i,2 = (F(i+1,1-F(i,1/(x0(i+2-x0(i; D(i,1 = 6 * F(i,2;
end %构建函数 D : d0 = 6*(F(1,2-f_0/h(1; ___________(4 dn = 6*(f_n-F(n-1,2/h(n-1; ___________(5 D = [d0;D;dn]; ______________(6 m= S\D; %寻找 x 所在位置,并求出对应插值: fori = 1:length(x for j = 1:n-1 if (x(i<=x0(j+1&(x(i>=x0(j y(i =( m(j*(x0(j+1-x(i^3/(6*h(j+... (m(j+1*(x(i-x0(j^3/(6*h(j+... (y0(j-(m(j*h(j^2/6*(x0(j+1-x(i/h(j+... (y0(j+1-(m(j+1*h(j^2/6*(x(i-x0(j/h(j ; break; else continue; end end end (2 (自然边界条件源代码: 仅仅需要对上面部分标注的位置做如下修改 :
(){} 2 1 ()(11),5,10,20: 1252 1()1,(0,1,2,,)()2,(0,1,2,,)() ()2 35,20:1100 (i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x = -≤≤=+=-+====-+ = 题目：插值多项式和三次样条插值多项式。已知对作、计算函数在点处的值；、求插值数据点 的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max ()n k n k n k n k n k n k k k N x S x k E N y N x E S y S x ==-=- 和； 、计算，； 解释你所得到的结果。 算法组织: 本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式 )(x N n 和三次样条插值多项式()n S x 。如此，则第三、四问则迎刃而解。计算两 种插值多项式的算法如下： 一、求Newton 插值多项式)(x N n ，算法组织如下： Newton 插值多项式的表达式如下： )())(()()(110010--???--+???+-+=n n n x x x x x x c x x c c x N 其中每一项的系数c i 的表达式如下： 1102110) ,,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -???-???= ???=- 根据i c 以上公式，计算的步骤如下： ?? ??? ?? ?????+??????? ???????????----) ,,,,(1) ,,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下：
题目一：多项式插值 某气象观测站在8：00（AM ）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton ）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。 二、数学原理 假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式： ）（） ）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα （1） 其中系数i α（i=0,1,2……n ）为特定系数，可由插值样条i i n y x P =） （（i=0,1,2……n ）确定。 根据均差的定义，把x 看成[a,b]上的一点,可得 f(x)= f （0x ）+f[10x x ，]（0x -x ） f[x, 0x ]= f[10x x ，]+f[x,10x x ，] （1x -x ） …… f[x, 0x ，…x 1-n ]= f[x, 0x ，…x n ]+ f[x, 0x ，…x n ](x-x n ) 综合以上式子，把后一式代入前一式，可得到： f(x)= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n )+ f[x, 0x ，…x n ,x ]） （x 1n +ω= N n （x ）+） （x n R 其中
N n （x ）= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n ) （2） ）（x n R = f(x)- N n （x ）= f[x, 0x ， (x) n ,x ]） （x 1n +ω （3） ） （x 1n +ω=（0x -x ）…(x-x n ) Newton 插值的系数i α（i=0,1,2……n ）可以用差商表示。一般有 f k =α[k 10x x x ??，] （k=0,1,2，……，n ） （4） 把（4）代入（1）得到满足插值条件N ）（） （i i n x f x =（i=0,1,2，……n ）的n 次Newton 插值多项式 N n （x ）=f （0x ）+f[10x x ，]（1x -x ）+f[210x x x ，，]（1x -x ）（2x -x ）+……+f[n 10x x x ??，]（1x -x ）（2x -x ）…（1-n x -x ）. 其中插值余项为： ） （）！ （） （）（）（）（x 1n f x N -x f x R 1n 1 n n +++==ωξ ξ介于k 10x x x ??，之间。 三、程序设计 function [y,A,C,L]=newdscg(X,Y,x,M) % y 为对应x 的值，A 为差商表，C 为多项式系数，L 为多项式 % X 为给定节点，Y 为节点值，x 为待求节点 n=length(X); m=length(x); % n 为X 的长度 for t=1:m
6.1 插值问题及其误差 6.1.2 与插值有关的MATLAB 函数 (一) POLY2SYM函数 调用格式一：poly2sym (C) 调用格式二：f1=poly2sym(C,'V') 或f2=poly2sym(C, sym ('V') ), (二) POLYVAL函数 调用格式：Y = polyval(P,X) (三) POLY函数 调用格式：Y = poly (V) (四) CONV函数 调用格式：C =conv (A, B) 例 6.1.2求三个一次多项式、和的积.它们的零点分别依次为0.4,0.8,1.2. 解我们可以用两种MATLAB程序求之. 方法1如输入MATLAB程序 >> X1=[0.4,0.8,1.2]; l1=poly(X1), L1=poly2sym (l1) 运行后输出结果为 l1 = 1.0000 - 2.4000 1.7600 -0.3840 L1 = x^3-12/5*x^2+44/25*x-48/125 方法2如输入MATLAB程序 >> P1=poly(0.4);P2=poly(0.8);P3=poly(1.2); C =conv (conv (P1, P2), P3) , L1=poly2sym (C) 运行后输出的结果与方法1相同. (五) DECONV 函数 调用格式：[Q,R] =deconv (B,A) (六) roots(poly(1:n))命令 调用格式：roots(poly(1:n)) (七) det(a*eye(size (A)) - A)命令 调用格式：b=det(a*ey e(size (A)) - A) 6.2 拉格朗日（Lagrange）插值及其MATLAB程序 6.2.1 线性插值及其MATLAB程序 例 6.2.1 已知函数在上具有二阶连续导数，，且满足条件 .求线性插值多项式和函数值，并估计其误差. 解输入程序 >> X=[1,3];Y=[1,2]; l01= poly(X(2))/( X(1)- X(2)), l11= poly(X(1))/( X(2)- X(1)), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x) 运行后输出基函数l0和l1及其插值多项式的系数向量P（略）、插值多项式L和插值Y为l0 = l1 = L = Y = -1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500 输入程序 >> M=5;R1=M*abs((x-X(1))* (x-X(2)))/2
ＭATLＡＢ三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（１，２，…,ｎ）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先，通过函数saｎwanｊ得到方程的系数矩阵，即追赶法方程的四个向量参数,接下来调 用追赶法(在inteｒｓaｎwj函数中），得到三次样条分段函数系数因子,然后进行多项式合并 得到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算 改点的值。附:追赶法程序cｈase ％%%％％%%%%%%%％% fｕncｔiｏｎ [nｅwv,w,newu，neｗd]=sanwj（ｘ,y,x０,y0,y1ａ,y1b）?％三弯矩样 条插值?％将插值点分两次输入，x0y０单独输入?% 边值条件a的二阶导数 y1a 和ｂ 的二阶导数ｙ1b，这里建议将y１ａ和y1b换成y2a和ｙ2b,以便于和三转角代码相区别 ?n=length(x);m=leｎgth(y); ｉf m~=n?errｏｒ(＇x or y 输入有误,再来'）; end?v=ｏneｓ(ｎ-１，1);u=ｏnes(n－１,1）;d=zeros（n－１,１）;?w=２*ｏ nes(n+１)；?h0=x（１）-x0;?h=zerｏs(n-１,1); for k＝1:n-１?ｈ(k)＝x（k+1)－x（k);?enｄ ｖ(1)＝h0/(ｈ0+h(1）); u(1)＝1－v（1）; ｄ（１)＝6*((y（2）-y(1）)/h（1)-(y(1)-y０)/ｈ0）/(h0+ｈ(１）);?% ｆoｒ k＝2:ｎ-1?ｖ(ｋ）＝h(ｋ-１)/(ｈ(k-１)+h(ｋ)）;?u(k)=1－v(k);?d（k)＝ 6*((y(ｋ+１)-ｙ（ｋ)）/h（k）-(ｙ(k）-y(ｋ-1）)/ｈ(k－1))／(h（k-1)+h(k)); ｅnd newv=[v;1］;?newu=［1;u]； d0=6＊（(y(1)-y0)/h0-y1a)/h0; d(n)=6＊（y1ｂ-(ｙ(n）－ｙ(n-1))／ｈ(n-1))／h(ｎ-1）; newd=［d0;d]; %%%%%%％％%％％％ ｆunctｉon iｎｔersａnｗj（x，y,ｘ0,y0,y1a，y1b) %三弯矩样条插值?%第一部分?n=leｎgtｈ(ｘ）;ｍ＝length(y); iｆ m~＝n?eｒror（'ｘｏｒy 输入有误，再来'); enｄ?%重新定义h?h=zｅｒｏs(n,1）; h(1)=x（1）-ｘ０; fｏr k=2:n h（k)=ｘ(ｋ)－x(k-1)；?end ％sptep1调用三弯矩函数?[a,b,c,d]=sanwj(x，ｙ,ｘ0,ｙ0,y１ａ,y1b);
(){} 21 ()(11),5,10,20: 1252 1()1,(0,1,2,,)()2,(0,1,2,,)() ()2 35,20:1100 (i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x =-≤≤=+=-+====-+ = 题目：插值多项式和三次样条插值多项式。 已知对作、计算函数在点处的值；、求插值数据点 的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max ()n k n k n k n k n k n k k k N x S x k E N y N x E S y S x ==-=- 和； 、计算，； 解释你所得到的结果。 算法组织: 本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式 )(x N n 和三次样条插值多项式()n S x 。如此，则第三、四问则迎刃而解。计算两种插值多项式的算法如下： 一、求Newton 插值多项式)(x N n ，算法组织如下： Newton 插值多项式的表达式如下： )())(()()(110010--???--+???+-+=n n n x x x x x x c x x c c x N 其中每一项的系数c i 的表达式如下： 1102110) ,,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -???-???= ???=- 根据i c 以上公式，计算的步骤如下： ?? ??? ?? ?????+??????? ???????????----) ,,,,(1) ,,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下：
数据插值和函数逼近 1数据插值 由已知样本点，以数据更为平滑为目标，求出其他点处的函数 值。在信号处理与图像处理上应用广泛。求解方法： yl=interpl (x, y, xl,'方法，) zl=in terp2 (x, y, z, xl, yl,'方法，) 1. 1 一维数据的插值 例：假设样本点来自/（兀）=（亍_3兀+ 50叫山兀,进行插值处理,得到平滑的曲线。 x二0:0. 1:1; y二(x?”2-3*x+5). *exp (-5*x). *sin(x); plot (x, y,' ; hold on ; xs二0:0. 01:1; yl= interpl (x, y, xs) ; % 线性插值 y2= interpl (x, y, xs, 'cubic'、 % 分段 3 次Hermite 扌西值y3= interpl (x, y, xs, 'spline'、% 分段三次样条插值plot (xs, yl,，r , xs, y2,，g , xs, y3,，b，); 例：草图样条曲线功能。 function sketch() x=[] ; y=[] ; gca; hold on; axis([0 1, 0 1]); while 1 [xO, yO, button]二ginput (1);
16s if (isempty (button)) break; end; x=[x xO]: y=[y yO]; plot (x, y,' ; end; xx 二[x(l) : (x(end)-x(1))/100:x(end)]; yy=interpl (x, y, xx,' spline'); plot (xx, yy, x, y ，‘*'); end 1.2二维网格数据的插值 例3 :假设样本点来自2 = (/_2兀0宀宀卩，进行插值处理，得到平滑的曲 面。 [x, y]=meshgrid(-3:0. 6:3, -2:0. 4:2); z 二(x.八 2 2*x)? * exp (-x ?"2—y. ^2~x ? *y) ; mesh (x, y, z); [xl, yl]=meshgrid(-3:0. 2:3, -2:0. 2:2); zl=interp2 (x, y, z, xl, yl) ;mesh(xl, yl, zl); zl=interp2 (x, y, z, xl, yl,' cubic') ;mesh (xl, yl, zl); zl=interp2 (x, y, z, xl, yl, ' spline') ;mesh(xl, yl, zl); / 4 -1 2 -2 -4 1.5 丫 1 . 0.5, 0. ?0.5 、 f ? K 、
告： Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi 是阶数为length(xi)*size(Y,2)的输出矩阵。 (2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函
数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数p chip，用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形； ’cubic’：与’pchip’操作相同； ’v5cubic’：在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。 (4)yi = interp1(x,Y,xi,method,'extrap') 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5)yi = interp1(x,Y,xi,method,extrapval) 确定超出x 范围的xi 中的分量的外插值extrapval，其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x); 2.>>xx = 0:.25:10; yy = interp1(x,y,xx); 3.>>plot(x,y,'kd',xx,yy) 复制代码 例2 1.>> year = 1900:10:2010; 2.>> product = [75.995 91.972 105.711 12 3.203 131.669 150.697 179.323 203.212 226.505