《matlab在科学计算中的应用062》由会员分享，可在线阅读，更多相关《matlab在科学计算中的应用062（52页珍藏版）》请在金锄头文库上搜索。
1、6.2.5样条插值的MATLAB表示样条函数的概念 样条(spline)：原指工程设计中使用的一种绘图工具(如富有 弹性的细木条或金属条),应用其将已知点连接成一条光滑曲线( 称为样条曲线),并使连接点处有连续的曲率.三次样条插值即由 此抽象出来. 样条函数：数学上将具有一定光滑性的分段多项式函数称 为样条函数. 具体的,给定区间a,b的一个划分:1 三次样条插值- 样条插值: 利用样条函数进行插值(例如分段线性插值是一次样条插值) - 三次样条插值: 已知函数f(x)在区间a,b上的n个节点的函数值, 求插值函数S(x)满足:2- 定义三次样条函数类： S=csapi(x,y)其中x=x1,
2、x2,.,xn, y=y1,y2,yn为样本点。S返回样条函数对象的插值结果，包括子区间点、 各区间点三次多项式系数等。- 可用 fnplt()绘制出插值结果，其调用格式： fnplt(S)- 对给定的向量xp, 可用fnval()函数计算, 其调用格式： yp=fnval(S,xp)其中得出的yp是xp上各点的插值结果。 Matlab 样条插值工具箱3 例：已知y0=sin(x0),x0= =0,0.4,1,2,pi, 求该函数三次样条插值结果 解法: x0=0,0.4,1,2,pi; y0=sin(x0); sp=csapi(x0,y0), fnplt(sp, r:); hold on,
3、sp = form: ppbreaks: 0 0.4000 1 2 3.1416coefs: 4x4 doublepieces: 4order: 4dim: 1 ezplot(sin(t),0,pi); plot(x0,y0,o)4- 在(0.4000, 1)区间内，插值多项式可以表示为 ：- 查看插值多项式 sp.coefs ans =-0.1627 0.0076 0.9965 0-0.1627 -0.1876 0.9245 0.38940.0244 -0.4804 0.5238 0.84150.0244 -0.4071 -0.3637 0.90935例: 点，用三次样条插值的方法对这些数据
4、进行拟合 解法: x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); sp=csapi(x,y); fnplt(sp)6Columns 7 through 12 0.4800 0.6000 -0.2404 0.7652 -0.5776 0.15880.6000 0.7200 -0.4774 0.6787 -0.4043 0.10010.7200 0.8400 -0.4559 0.5068 -0.2621 0.06050.8400 0.9600 -0.4559 0.3427 -0.1601 0.0356 c=sp.breaks(1:4) sp.breaks(
5、2:5) sp.coefs(1:4,:),. sp.breaks(5:8) sp.breaks(6:9) sp.coefs(5:8,:) c =Columns 1 through 6 0 0.1200 24.7396 -19.3588 4.5151 0 0.1200 0.2400 24.7396 -10.4526 0.9377 0.3058 0.2400 0.3600 4.5071 -1.5463 -0.5022 0.31050.3600 0.4800 1.9139 0.0762 -0.6786 0.23587- 格式S=csapi(x1,x2,xn,z)注: csapi() 可处理多个自变量
6、的网格数据三次样条插值类8 x0=-3:.6:3; y0=-2:.4:2; x,y=ndgrid(x0,y0); % 注意这里 只能用 ndgrid, 否则生成的 z 矩阵顺序有问题 z=(x.2-2*x).* exp(-x.2-y.2-x.*y); sp=csapi(x0,y0,z); fnplt(sp);例:对函数解法:9 三次样条数据插值 -Matlab实现: 函数spline 格式 yy = spline(x,y,xx) 例:对离散分布在y=exp(x)sin(x)函数曲线上的 数据点进行样条插值计算：x = 0 2 4 5 8 12 12.8 17.2 19.9 20; y = ex
7、p(x).*sin(x);xx = 0:.25:20;yy = spline(x,y,xx);plot(x,y,o,xx,yy)1011B样条函数及其MATLAB表示S=spapi(k,x,y)- Matlab实现:12 x0=0,0.4,1,2,pi; y0=sin(x0); ezplot(sin(t),0,pi); hold on sp1=csapi(x0,y0); fnplt(sp1,-); % 三次分段多项式样条插值例:y=sin(t)和解法: sp2=spapi(5,x0,y0); fnplt(sp2,r:) % 5 次 B 样条插值13 x=0:.12:1; y=(x.2-3*x+
8、5).*exp(-5*x).*sin(x); ezplot(x2-3*x+5)*exp(-5*x)*sin(x),0,1), hold on sp1=csapi(x,y); fnplt(sp1,-); sp2=spapi(5,x,y);fnplt(sp2,r:)B样条的插值效果好于三次分段多项式146.2.6 基于样条插值的数值微积分运算 基于样条插值的数值微分运算 格式：Sd=fnder(S,k) 该函数可以求取S的k阶导数。格式：Sd=fnder(S,k1,kn) 可以求取多变量函数的偏导数15 例：解法： syms x; f=(x2-3*x+5)*exp(-5*x)*sin(x); ez
9、plot(diff(f),0,1), hold on 解析解微分图 x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x);取样本点 sp1=csapi(x,y);建立三次样条函数 dsp1=fnder(sp1,1);数值微分 fnplt(dsp1, r-)绘制样条图 sp2=spapi(5,x,y);5阶次B样条 dsp2=fnder(sp2,1); 数值微分 fnplt(dsp2,g:); axis(0,1,-0.8,5)1617 例：解法： Step1：拟合曲面 x0=-3:.3:3; y0=-2:.2:2; x,y=ndgrid(x0,y0); z=(x
10、.2-2*x).*exp(-x.2-y.2-x.*y); sp=spapi(5,5, x0,y0,z); B样条 dspxy=fnder(sp,1,1); fnplt(dspxy) %生成样条图18理论方法： syms x y; z=(x2-2*x)*exp(-x2-y2-x*y); ezsurf(diff(diff(z,x),y),-3 3,-2 2)对符号变量表达式做三维表面图19 基于样条插值的数值积分运算 格式：f=fnint(S) 其中S为样条函数。解法： x=0,0.4,1 2,pi; y=sin(x); sp1=csapi(x,y); a=fnint(sp1); 建立三次样条函数
11、并积分 xx=fnval(a,0,pi); xx(2)-xx(1) ans =2.0191 例：考虑中较稀疏的样本点,用样条积分的方式求出定积分及积分函数。20Step2: sp2=spapi(5,x,y); b=fnint(sp2); 建立B样条函数并积分 xx=fnval(b,0,pi); xx(2)-xx(1) ans =1.9999 Step3: 绘制曲线 ezplot(-cos(t)+2,0,pi);hold on fnplt(a,r-); fnplt(b,g:)216.3 数据拟合 用插值的方法对一函数进行近似,要求所得到的 插值多项式经过已知插值节点;在n比较大的情况 下,插值多
12、项式往往是高次多项式,这也就容易出 现振荡现象（Runge现象），即虽然在插值节点 上没有误差,但在插值节点之外插值误差变得很 大,从“整体”上看,插值逼近效果将变得“很差”。数据拟合, 又称函数逼近, 求近似函数的一类数值方法。 所谓数据拟合是求一个简单的函数,例如是一个 低次多项式,不要求通过已知的这些点,而是要求 在整体上“尽量好”的逼近原函数。这时,在每个 已知点上就会有误差,数据拟合就是从整体上使 误差,尽量的小一些。22 例：给定一组实验数据,如下表 i1432 2648 1.12.84.97.2086422468求x，y的函数关系。解：做草图(见左图), 由此可设困难: 无论a1
13、,a0取何值,直线都不 能经过所有样本点.问题: 选取适当a1,a0的值,使直线能最 好的反映数据点的基本趋势.23 残差：近似值与样本点函数值（已知数据）之差 常用准则：1). 残差绝对值之和最小：2). 残差最大绝对值最小：3). 残差平方和最小：1)使用不方便； 2) 准则求近似函数的方法称为函数的最佳一致逼近 3) 准则求近似函数的方法称为最佳平方逼近，或数据拟合的 最小二乘法246.3.1 多项式拟合 n次多项式： 已知数据点： 曲线与数据点的残差为： 最小二乘法函数逼近，即取适当 使得最小.由极值原理，得方程组： 数学算法25 或 或矩阵形式：26多项式拟合MATLAB命令：pol
14、yfit格式：p=polyfit(x,y,n) Matlab实现：27Step1:多项式拟合 x0=0:.1:1; y0=(x0.2-3*x0+5).*exp(- 5*x0).*sin(x0); p3=polyfit(x0,y0,3); vpa(poly2sym(p3),10) % 可以如下显示多项式 ans = 2.839962923*x3- 4.789842696*x2+1.943211631*x+.5975248921e-1例：解法：28Step2：绘制拟合曲线并与精确解曲线比较： x=0:.01:1; ya=(x.2-3*x+5).*exp(-5*x).*sin(x); y1=polyval(p3,x); plot(x,y1,x,ya,x0,y0,o) legend(拟合曲线,精确解曲线,样本点)29Step3：就不同的次数进行拟合： p4=polyfit(x0,y0,4); y2=polyval(p4,x); p5=polyfit(x0,y0,5); y3=polyval(p5,x); p8=polyfit(x0,y0,8); y4=polyval(p8,x); plot(x,ya,x0,y0,o,x,y2,x,y3,x,y4)最小二乘8次 多项式拟合 效
