matlab在科学计算中的应用.ppt

1,6.2.5样条插值的MATLAB表示,样条函数的概念,样条(spline)：原指工程设计中使用的一种绘图工具(如富有弹性的细木条或金属条),应用其将已知点连接成一条光滑曲线(称为样条曲线),并使连接点处有连续的曲率.三次样条插值即由此抽象出来.,样条函数：数学上将具有一定光滑性的分段多项式函数称为样条函数.具体的,给定区间a,b的一个划分:,2,三次样条插值-样条插值:利用样条函数进行插值(例如分段线性插值是一次样条插值),-三次样条插值:已知函数f(x)在区间a,b上的n个节点,3,-定义三次样条函数类：S=csapi(x,y)其中x=x1,x2,.,xn,y=y1,y2,yn为样本点。S返回样条函数对象的插值结果，包括子区间点、各区间点三次多项式系数等。-可用fnplt()绘制出插值结果，其调用格式：fnplt(S)-对给定的向量xp,可用fnval()函数计算,其调用格式：yp=fnval(S,xp)其中得出的yp是xp上各点的插值结果。,Matlab样条插值工具箱,4,例：已知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:);holdon,sp=form:ppbreaks:00.4000123.1416coefs:4x4doublepieces:4order:4dim:1ezplot(sin(t),0,pi);plot(x0,y0,o),5,-在(0.4000,1)区间内，插值多项式可以表示为：,-查看插值多项式sp.coefsans=-0.16270.00760.99650-0.1627-0.18760.92450.38940.0244-0.48040.52380.84150.0244-0.4071-0.36370.9093,6,例:,点，用三次样条插值的方法对这些数据进行拟合,解法:x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);sp=csapi(x,y);fnplt(sp),7,Columns7through120.48000.6000-0.24040.7652-0.57760.15880.60000.7200-0.47740.6787-0.40430.10010.72000.8400-0.45590.5068-0.26210.06050.84000.9600-0.45590.3427-0.16010.0356,c=sp.breaks(1:4)sp.breaks(2:5)sp.coefs(1:4,:),.sp.breaks(5:8)sp.breaks(6:9)sp.coefs(5:8,:)c=Columns1through600.120024.7396-19.35884.515100.12000.240024.7396-10.45260.93770.30580.24000.36004.5071-1.5463-0.50220.31050.36000.48001.91390.0762-0.67860.2358,8,-格式S=csapi(x1,x2,xn,z),注:csapi()可处理多个自变量的网格数据三次样条插值类,9,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);,例:对函数,解法:,10,三次样条数据插值-Matlab实现:函数spline格式yy=spline(x,y,xx)例:对离散分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算：x=024581212.817.219.920;y=exp(x).*sin(x);xx=0:.25:20;yy=spline(x,y,xx);plot(x,y,o,xx,yy),11,12,B样条函数及其MATLAB表示,S=spapi(k,x,y),-Matlab实现:,13,x0=0,0.4,1,2,pi;y0=sin(x0);ezplot(sin(t),0,pi);holdonsp1=csapi(x0,y0);fnplt(sp1,-);%三次分段多项式样条插值,解法:,sp2=spapi(5,x0,y0);fnplt(sp2,r:)%5次B样条插值,14,x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);ezplot(x2-3*x+5)*exp(-5*x)*sin(x),0,1),holdonsp1=csapi(x,y);fnplt(sp1,-);sp2=spapi(5,x,y);fnplt(sp2,r:),B样条的插值效果好于三次分段多项式,15,6.2.6基于样条插值的数值微积分运算,基于样条插值的数值微分运算格式：Sd=fnder(S,k)该函数可以求取S的k阶导数。格式：Sd=fnder(S,k1,kn)可以求取多变量函数的偏导数,16,例：,解法：symsx;f=(x2-3*x+5)*exp(-5*x)*sin(x);ezplot(diff(f),0,1),holdon解析解微分图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),17,18,例：解法：Step1：拟合曲面x0=-3:.3:3;y0=-2:.2:2;x,y=ndgrid(x0,y0);z=(x.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)%生成样条图,19,理论方法：symsxy;z=(x2-2*x)*exp(-x2-y2-x*y);ezsurf(diff(diff(z,x),y),-33,-22)对符号变量表达式做三维表面图,20,基于样条插值的数值积分运算格式：f=fnint(S)其中S为样条函数。,解法：x=0,0.4,12,pi;y=sin(x);sp1=csapi(x,y);a=fnint(sp1);建立三次样条函数并积分xx=fnval(a,0,pi);xx(2)-xx(1)ans=2.0191,用样条积分的方式求出定积分及积分函数。,21,Step2:sp2=spapi(5,x,y);b=fnint(sp2);建立B样条函数并积分xx=fnval(b,0,pi);xx(2)-xx(1)ans=1.9999Step3:绘制曲线ezplot(-cos(t)+2,0,pi);holdonfnplt(a,r-);fnplt(b,g:),22,6.3数据拟合,用插值的方法对一函数进行近似,要求所得到的插值多项式经过已知插值节点;在n比较大的情况下,插值多项式往往是高次多项式,这也就容易出现振荡现象（Runge现象），即虽然在插值节点上没有误差,但在插值节点之外插值误差变得很大,从“整体”上看,插值逼近效果将变得“很差”。,数据拟合,又称函数逼近,求近似函数的一类数值方法。,所谓数据拟合是求一个简单的函数,例如是一个低次多项式,不要求通过已知的这些点,而是要求在整体上“尽量好”的逼近原函数。这时,在每个已知点上就会有误差,数据拟合就是从整体上使误差,尽量的小一些。,23,例：给定一组实验数据,如下表,求x，y的函数关系。,解：做草图(见左图),由此可设,困难:无论a1,a0取何值,直线都不能经过所有样本点.,问题:选取适当a1,a0的值,使直线能最好的反映数据点的基本趋势.,24,残差：近似值与样本点函数值（已知数据）之差,常用准则：,1).残差绝对值之和最小：,2).残差最大绝对值最小：,3).残差平方和最小：,1)使用不方便；2)准则求近似函数的方法称为函数的最佳一致逼近3)准则求近似函数的方法称为最佳平方逼近，或数据拟合的最小二乘法,25,6.3.1多项式拟合,n次多项式：,已知数据点：,曲线与数据点的残差为：,最小二乘法函数逼近，即取适当使得,最小.由极值原理，得方程组：,数学算法,26,或,或矩阵形式：,27,多项式拟合MATLAB命令：polyfit格式：p=polyfit(x,y,n),Matlab实现：,28,Step1:多项式拟合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,例：,解法：,29,Step2：绘制拟合曲线并与精确解曲线比较：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(拟合曲线,精确解曲线,样本点),30,Step3：就不同的次数进行拟合：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次多项式拟合效果较好,31,注：察看拟合最高次数为8的多项式：vpa(poly2sym(p8),5)ans=-8.2586*x8+43.566*x7-101.98*x6+140.22*x5-125.29*x4+74.450*x3-27.672*x2+4.9869*x+.42037e-6Taylor幂级数展开：symsx;y=(x2-3*x+5)*exp(-5*x)*sin(x);vpa(taylor(y,9),5)ans=5.*x-28.*x2+77.667*x3-142.*x4+192.17*x5-204.96*x6+179.13*x7-131.67*x8,多项式表示数据模型是不唯一的，即是两个多项式函数完全不同。在某一区域内其曲线可能特别近似,32,解法：Step1：取样本点和精确数值解x0=-1+2*0:10/10;y0=1./(1+25*x0.2);x=-1:.01:1;ya=1./(1+25*x.2);Step2：各次多项式拟合p3=polyfit(x0,y0,3);y1=polyval(p3,x);p5=polyfit(x0,y0,5);y2=polyval(p5,x);p8=polyfit(x0,y0,8);y3=polyval(p8,x);p10=polyfit(x0,y0,10);,例：,33,y4=polyval(p10,x);plot(x,ya,x,y1,x,y2,-.,x,y3,-,x,y4,:)legend(精确解曲线,3次多项式拟合,5次多项式拟合,8次多项式拟合,10次多项式拟合),多项式拟合的效果并不一定总是很精确的。,34,考察Taylor幂级数展开拟合效果symsx;y=1/(1+25*x2);p=taylor(y,x,10)p=1-25*x2+625*x4-15625*x6+390625*x8x1=-1:0.01:1;ya=1./(1+25*x1.2);y1=subs(p,x,x1);plot(x1,ya,-,x1,y1),Taylor多项式拟合效果更差,35,6.3.2函数线性组合的曲线拟合方法,且线性无关,36,该方程的最小二乘解为:,其中,37,例,38,x=0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5;y=2.88;2.2576;1.9683;1.9258;2.0862;2.109;2.1979;2.5409;2.9627;3.155;3.2052;A=ones(size(x)，exp(-3*x),cos(-2*x).*exp(-4*x)，x.2;c=Ay;c1=cc1=1.22002.3397-0.67970.8700,39,图形显示x0=0:0.01:1.5;A1=ones(size(x0)exp(-3*x0),cos(-2*x0).*exp(-4*x0)x0.2;y1=A1*c;plot(x0,y1,x,y,x),40,数据分析x=1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,.2.2255,2.4596,2.7183,3.6693;y=0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,.0.2864,0.2532,0.2238,0.1546;plot(x,y,x,y,*),例,41,分别对x,y进行对数变换：x1=log(x);y1=log(y);plot(x1,y1),数据点xi和yi的分布近似于指数曲线，那么数组(xi，ln(yi)的分布应近似于直线,42,A=x1,ones(size(x1);c=Ay1c=-1.2339-0.2630exp(c(2)ans=0.7687,43,解法：x=0:0.1:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);n=8;A=;fori=1:n+1,A(:,i)=x.(n+1-i);endc=Ay;vpa(poly2sym(c),5)ans=-8.2586*x8+43.566*x7-101.98*x6+140.22*x5-125.29*x4+74.450*x3-27.672*x2+4.9869*x+.42037e-6,例,与polyfit作用结果一致,44,6.3.3最小二乘曲线拟合,数学描述：,45,格式：a,jm=lsqcurvefit(Fun,a0,x,y),Matlab实现：,46,例：,f=inline(a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x),a,x);,解法：,Step1:,47,xx,res=lsqcurvefit(f,1,1,1,1,1,x,y);xx,resOptimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOP