《Matlab数据处理.ppt》由会员分享，可在线阅读，更多相关《Matlab数据处理.ppt（131页珍藏版）》请在三一办公上搜索。
1、第四章 MATLAB数据处理,声明和致谢,课件取自网络电子课件中科院MATLAB课件，只进行了很少的编辑和加工，对不知名的网络电子课件作者表示感谢！,主要内容,微积分问题的解析解函数的级数展开与级数求和问题求解概率分布与伪随机数生成数据插值数据拟合,4.1 微积分问题的解析解 4.1.1 极限问题的解析解,单变量函数的极限格式1：L=limit(fun,x,x0)格式2：L=limit(fun,x,x0,left 或 right),例:试求解极限问题 syms x a b;f=x*(1+a/x)x*sin(b/x);L=limit(f,x,inf)L=exp(a)*b例:求解单边极限问题 sy
2、ms x;limit(exp(x3)-1)/(1-cos(sqrt(x-sin(x),x,0,right)ans=12,在(-0.1,0.1)区间绘制出函数曲线：x=-0.1:0.001:0.1;y=(exp(x.3)-1)./(1-cos(sqrt(x-sin(x);Warning:Divide by zero.(Type warning off MATLAB:divideByZero to suppress this warning.)plot(x,y,-,0,12,o),4.1.2 函数导数的解析解,函数的导数和高阶导数格式：y=diff(fun,x)%求导数(默认为1阶)y=diff(
3、fun,x,n)%求n阶导数例：一阶导数：syms x;f=sin(x)/(x2+4*x+3);f1=diff(f);pretty(f1),cos(x)sin(x)(2 x+4)-2 2 2 x+4 x+3(x+4 x+3)原函数及一阶导数图：x1=0:.01:5;y=subs(f,x,x1);y1=subs(f1,x,x1);plot(x1,y,x1,y1,:)更高阶导数：tic,diff(f,x,100);tocelapsed_time=4.6860,原函数4阶导数 f4=diff(f,x,4);pretty(f4)2 sin(x)cos(x)(2 x+4)sin(x)(2 x+4)-+4
4、-12-2 2 2 2 3 x+4 x+3(x+4 x+3)(x+4 x+3)3 sin(x)cos(x)(2 x+4)cos(x)(2 x+4)+12-24-+48-2 2 2 4 2 3(x+4 x+3)(x+4 x+3)(x+4 x+3)4 2 sin(x)(2 x+4)sin(x)(2 x+4)sin(x)+24-72-+24-2 5 2 4 2 3(x+4 x+3)(x+4 x+3)(x+4 x+3),多元函数的偏导：格式：f=diff(diff(f,x,m),y,n)或 f=diff(diff(f,y,n),x,m)例：求其偏导数并用图表示。syms x y z=(x2-2*x)*
5、exp(-x2-y2-x*y);zx=simple(diff(z,x)zx=-exp(-x2-y2-x*y)*(-2*x+2+2*x3+x2*y-4*x2-2*x*y),zy=diff(z,y)zy=(x2-2*x)*(-2*y-x)*exp(-x2-y2-x*y)直接绘制三维曲面 x,y=meshgrid(-3:.2:3,-2:.2:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);surf(x,y,z),axis(-3 3-2 2-0.7 1.5),contour(x,y,z,30),hold on%绘制等值线 zx=-exp(-x.2-y.2-x.*y).*(-2*x
6、+2+2*x.3+x.2.*y-4*x.2-2*x.*y);zy=-x.*(x-2).*(2*y+x).*exp(-x.2-y.2-x.*y);%偏导的数值解 quiver(x,y,zx,zy)%绘制引力线,例 syms x y z;f=sin(x2*y)*exp(-x2*y-z2);df=diff(diff(diff(f,x,2),y),z);df=simple(df);pretty(df)2 2 2 2 2-4 z exp(-x y-z)(cos(x y)-10 cos(x y)y x+4 2 4 2 2 4 2 2sin(x y)x y+4 cos(x y)x y-sin(x y),4.
7、1.3 积分问题的解析解,不定积分的推导：格式：F=int(fun,x)例：用diff()函数求其一阶导数，再积分，检验是否可以得出一致的结果。syms x;y=sin(x)/(x2+4*x+3);y1=diff(y);y0=int(y1);pretty(y0)%对导数积分 sin(x)sin(x)-1/2-+1/2-x+3 x+1,对原函数求4 阶导数，再对结果进行4次积分 y4=diff(y,4);y0=int(int(int(int(y4);pretty(simple(y0)sin(x)-2 x+4 x+3,例:证明 syms a x;f=simple(int(x3*cos(a*x)2,
8、x)f=1/16*(4*a3*x3*sin(2*a*x)+2*a4*x4+6*a2*x2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a4 f1=x4/8+(x3/(4*a)-3*x/(8*a3)*sin(2*a*x)+.(3*x2/(8*a2)-3/(16*a4)*cos(2*a*x);simple(f-f1)%求两个结果的差ans=-3/16/a4,定积分与无穷积分计算:格式：I=int(f,x,a,b)格式：I=int(f,x,a,inf),例：syms x;I1=int(exp(-x2/2),x,0,1.5)无解I1=1/2*erf(3/4*
9、2(1/2)*2(1/2)*pi(1/2)vpa(I1,70)ans=I2=int(exp(-x2/2),x,0,inf)I2=1/2*2(1/2)*pi(1/2),多重积分问题的MATLAB求解例：syms x y z;f0=-4*z*exp(-x2*y-z2)*(cos(x2*y)-10*cos(x2*y)*y*x2+.4*sin(x2*y)*x4*y2+4*cos(x2*y)*x4*y2-sin(x2*y);f1=int(f0,z);f1=int(f1,y);f1=int(f1,x);f1=simple(int(f1,x)f1=exp(-x2*y-z2)*sin(x2*y),f2=int
10、(f0,z);f2=int(f2,x);f2=int(f2,x);f2=simple(int(f2,y)f2=2*exp(-x2*y-z2)*tan(1/2*x2*y)/(1+tan(1/2*x2*y)2)simple(f1-f2)ans=0 顺序的改变使化简结果不同于原函数，但其误差为0，表明二者实际完全一致。这是由于积分顺序不同，得不出实际的最简形式。,例：syms x y z int(int(int(4*x*z*exp(-x2*y-z2),x,0,1),y,0,pi),z,0,pi)ans=(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2)*pi2*hyper
11、geom(1,2,-pi2)Ei(n,z)为指数积分，无解析解，但可求其数值解：vpa(ans,60)ans=,主要内容,微积分问题的解析解函数的级数展开与级数求和问题求解概率分布与伪随机数生成数据插值数据拟合,4.2 函数的级数展开与 级数求和问题求解,4.2.1 Taylor 幂级数展开4.2.2 Fourier 级数展开4.2.3 级数求和的计算,4.2.1 Taylor 幂级数展开 4.2.1.1 单变量函数的 Taylor 幂级数展开,例:syms x;f=sin(x)/(x2+4*x+3);y1=taylor(f,x,9);pretty(y1)2 23 3 34 4 4087 5
12、3067 6 515273 7 386459 8 1/3 x-4/9 x+-x-x+-x-x+-x-x 54 81 9720 7290 1224720 918540,taylor(f,x,9,2)ans=syms a;taylor(f,x,5,a)%结果较冗长，显示从略ans=sin(a)/(a2+3+4*a)+(cos(a)-sin(a)/(a2+3+4*a)*(4+2*a)/(a2+3+4*a)*(x-a)+(-sin(a)/(a2+3+4*a)-1/2*sin(a)-(cos(a)*a2+3*cos(a)+4*cos(a)*a-4*sin(a)-2*sin(a)*a)/(a2+3+4*a
13、)2*(4+2*a)/(a2+3+4*a)*(x-a)2+,例:对y=sinx进行Taylor幂级数展开,并观察不同阶次的近似效果。x0=-2*pi:0.01:2*pi;y0=sin(x0);syms x;y=sin(x);plot(x0,y0,r-.),axis(-2*pi,2*pi,-1.5,1.5);hold on for n=8:2:16 p=taylor(y,x,n),y1=subs(p,x,x0);line(x0,y1);endp=x-1/6*x3+1/120*x5-1/5040*x7p=x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9p=x-1/6*
14、x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11p=x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11+1/6227020800*x13,p=,4.2.2 级数求和的计算,是在符号工具箱中提供的,例:计算 format long;sum(2.0:63)%数值计算ans=1.844674407370955e+019 sum(sym(2).0:200)%或 syms k;symsum(2k,0,200)把2定义为符号量可使计算更精确ans=syms k;symsum(2k,0,200)ans=
15、,例:试求解无穷级数的和 syms n;s=symsum(1/(3*n-2)*(3*n+1),n,1,inf)%采用符号运算工具箱s=1/3 m=1:10000000;s1=sum(1./(3*m-2).*(3*m+1);%数值计算方法,双精度有效位16,“大数吃小数”,无法精确 format long;s1%以长型方式显示得出的结果s1=0.33333332222165,例:求解 syms n x s1=symsum(2/(2*n+1)*(2*x+1)(2*n+1),n,0,inf);simple(s1)%对结果进行化简，MATLAB 6.5 及以前版本因本身 bug 化简很麻烦ans=lo
16、g(2*x+1)2)(1/2)+1)/(2*x+1)2)(1/2)-1)%实际应为log(x+1)/x),例:求 syms m n;limit(symsum(1/m,m,1,n)-log(n),n,inf)ans=eulergamma vpa(ans,70)%显示 70 位有效数字ans=,4.2.3 二元函数的梯度计算,格式：若z矩阵是建立在等间距的形式生成的网格基础上，则实际梯度为,例：计算梯度，绘制引力线图：x,y=meshgrid(-3:.2:3,-2:.2:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);fx,fy=gradient(z);fx=fx/0.2;f
17、y=fy/0.2;contour(x,y,z,50);hold on;quiver(x,y,fx,fy)%绘制等高线与引力线图,绘制误差曲面:zx=-exp(-x.2-y.2-x.*y).*(-2*x+2+2*x.3+x.2.*y-4*x.2-2*x.*y);zy=-x.*(x-2).*(2*y+x).*exp(-x.2-y.2-x.*y);surf(x,y,abs(fx-zx);axis(-3 3-2 2 0,0.08)figure;surf(x,y,abs(fy-zy);axis(-3 3-2 2 0,0.11)建立一个新图形窗口,为减少误差，对网格加密一倍：x,y=meshgrid(-3
18、:.1:3,-2:.1:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);fx,fy=gradient(z);fx=fx/0.1;fy=fy/0.1;zx=-exp(-x.2-y.2-x.*y).*(-2*x+2+2*x.3+x.2.*y-4*x.2-2*x.*y);zy=-x.*(x-2).*(2*y+x).*exp(-x.2-y.2-x.*y);surf(x,y,abs(fx-zx);axis(-3 3-2 2 0,0.02)figure;surf(x,y,abs(fy-zy);axis(-3 3-2 2 0,0.06),主要内容,微积分问题的解析解函数的级数展开与级数
19、求和问题求解概率分布与伪随机数生成数据插值数据拟合,4.3概率分布与伪随机数生成 4.3.1 概率密度函数与分布函数概述,通用函数计算概率密度函数值,函数 pdf格式 P=pdf(name，K，A)P=pdf(name，K，A，B)P=pdf(name，K，A，B，C)说明 返回在X=K处、参数为A、B、C的概率密度值，对于不同的分布，参数个数是不同；name为分布函数名。例如二项分布：设一次试验，事件Y发生的概率为p，那么，在n次独立重复试验中，事件Y恰好发生K次的概率P_K为：P_K=PX=K=pdf(bino，K，n，p),例:计算正态分布N（0，1）的随机变量X在点0.6578的密度函
20、数值。解：pdf(norm,0.6578,0,1)ans=0.3213例:自由度为8的卡方分布，在点2.18处的密度函数值。解：pdf(chi2,2.18,8)ans=0.0363,随机变量的累积概率值(分布函数值),通用函数cdf用来计算随机变量的概率之和（累积概率值）函数 cdf格式 cdf(name，K，A)cdf(name，K，A，B)cdf(name，K，A，B，C)说明 返回以name为分布、随机变量XK的概率之和的累积概率值，name为分布函数名.,例:求标准正态分布随机变量X落在区间(-，0.4)内的概率。解：cdf(norm,0.4,0,1)ans=0.6554例：求自由度为
21、16的卡方分布随机变量落在0，6.91内的概率。解：cdf(chi2,6.91,16)ans=0.0250,随机变量的逆累积分布函数,MATLAB中的逆累积分布函数是已知，求x。命令 icdf 计算逆累积分布函数格式 icdf(name，K，A)icdf(name，K，A，B)icdf(name，K，A，B，C)说明 返回分布为name，参数为a1,a2,a3，累积概率值为P的临界值，这里name与前面相同。如果F=cdf(name，X，A，B，C)，则 X=icdf(name，F，A，B，C),例:在标准正态分布表中，若已知F=0.6554,求X解：icdf(norm,0.6554,0,1)
22、ans=0.3999例：公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X（单位：cm）服从正态分布N（175，6），求车门的最低高度。解：设h为车门高度，X为身高。求满足条件 FXh=0.01故 h=icdf(norm,0.99,175,6)h=188.9581,4.3.2 常见分布的概率密度函数与分布函数 4.3.2.1 Poisson分布,其要求x是正整数。,其中：x为选定的一组横坐标向量，y为x各点处的概率密度函数值。,例：绘制 l=1,2,5,10 时 Poisson 分布的概率密度函数与概率分布函数曲线。x=0:15;y1=;y2=;lam1=1,2,5,
23、10;for i=1:length(lam1)y1=y1,poisspdf(x,lam1(i);y2=y2,poisscdf(x,lam1(i);end plot(x,y1),figure;plot(x,y2),4.3.2.2 正态分布,正态分布的概率密度函数为:,例：x=-5:.02:5;y1=;y2=;mu1=-1,0,0,0,1;sig1=1,0.1,1,10,1;sig1=sqrt(sig1);for i=1:length(mu1)y1=y1,normpdf(x,mu1(i),sig1(i);y2=y2,normcdf(x,mu1(i),sig1(i);end plot(x,y1),f
24、igure;plot(x,y2),4.3.2.3 分布,例：x=-0.5:.02:5;x=-eps:-0.02:-0.5,0:0.02:5;x=sort(x);替代 y1=;y2=;a1=1,1,2,1,3;lam1=1,0.5,1,2,1;for i=1:length(a1)y1=y1,gampdf(x,a1(i),lam1(i);y2=y2,gamcdf(x,a1(i),lam1(i);end plot(x,y1),figure;plot(x,y2),4.3.2.4 分布（卡方分布）,其为一特殊的 分布，a=k/2,l=1/2。,例：x=-eps:-0.02:-0.5,0:0.02:2;x
25、=sort(x);k1=1,2,3,4,5;y1=;y2=;for i=1:length(k1)y1=y1,chi2pdf(x,k1(i);y2=y2,chi2cdf(x,k1(i);end plot(x,y1),figure;plot(x,y2),4.3.2.5 分布,概率密度函数为:,其为参数k的函数，且k为正整数。,例：x=-5:0.02:5;k1=1,2,5,10;y1=;y2=;for i=1:length(k1)y1=y1,tpdf(x,k1(i);y2=y2,tcdf(x,k1(i);end plot(x,y1),figure;plot(x,y2),4.3.2.6 Rayleig
26、h分布,例：x=-eps:-0.02:-0.5,0:0.02:5;x=sort(x);b1=.5,1,3,5;y1=;y2=;for i=1:length(b1)y1=y1,raylpdf(x,b1(i);y2=y2,raylcdf(x,b1(i);end plot(x,y1),figure;plot(x,y2),4.3.2.7 F 分布,其为参数p,q的函数，且p,q均为正整数。,例：分别绘制（p,q）为（1,1),(2,1),(3,1)(3,2),(4,1)时F分布的概率密度函数与分布函数曲线。x=-eps:-0.02:-0.5,0:0.02:1;x=sort(x);p1=1 2 3 3
27、4;q1=1 1 1 2 1;y1=;y2=;for i=1:length(p1)y1=y1,fpdf(x,p1(i),q1(i);y2=y2,fcdf(x,p1(i),q1(i);end plot(x,y1),figure;plot(x,y2),4.3.3 概率问题的求解,图4-9,例：b=1;p1=raylcdf(0.2,b);p2=raylcdf(2,b);P1=p2-p1P1=0.8449 p1=raylcdf(1,b);P2=1-p1P2=0.6065,例：syms x y;f=x2+x*y/3;P=int(int(f,x,0,1/2),y,0,1/2)P=5/192 syms x
28、y;f=x2+x*y/3;P=int(int(f,x,0,1),y,0,2)P=1,4.3.4 随机数与伪随机数,例：b=1;p=raylrnd(1,30000,1);xx=0:.1:4;yy=hist(p,xx);hist()找出随机数落入各个子区间的点个数，并由之拟合出生成数据的概率密度。yy=yy/(30000*0.1);bar(xx,yy),y=raylpdf(xx,1);line(xx,y),主要内容,微积分问题的解析解函数的级数展开与级数求和问题求解概率分布与伪随机数生成数据插值数据拟合,一个多项式的幂级数形式可表示为：也可表为嵌套形式或因子形式 N阶多项式n个根，其中包含重根和复
29、根。若多项式所有系数均为实数，则全部复根都将以共轭对的形式出现,4.4 数据插值关于多项式MATLAB命令,幂系数：在MATLAB里，多项式用行向量表示，其元素为多项式的系数，并从左至右按降幂排列。例：被表示为 p=2 1 4 5 poly2sym(p)ans=2*x3+x2+4*x+5Roots:多项式的零点可用命令roots求的。例：r=roots(p)得到 r=0.2500+1.5612i 0.2500-1.5612i-1.0000 所有零点由一个列向量给出。,poly:由零点可得原始多项式的各系数，但可能相差一个常数倍。例：poly(r)ans=1.0000 0.5000 2.0000
30、 2.5000 注意：若存在重根，这种转换可能会降低精度。例：r=roots(1-6 15-20 15-6 1)r=1.0042+0.0025i 1.0042-0.0025i 1.0000+0.0049i 1.0000-0.0049i 0.9958+0.0024i 0.9958-0.0024i舍入误差的影响，与计算精度有关。,polyval:可用命令polyval计算多项式的值。例：计算y(2.5)c=3,-7,2,1,1；xi=2.5;yi=polyval(c,xi)yi=23.8125如果xi是含有多个横坐标值的数组，则yi也为与xi长度相同的向量。c=3,-7,2,1,1;xi=2.5,
31、3;yi=polyval(c,xi)yi=23.8125 76.0000,polyfit:给定n+1个点将可以唯一确定一个n阶多项式。利用命令polyfit可容易确定多项式的系数。例：x=1.1,2.3,3.9,5.1;y=3.887,4.276,4.651,2.117;a=polyfit(x,y,length(x)-1)a=-0.2015 1.4385-2.7477 5.4370 poly2sym(a)ans=-403/2000*x3+2877/2000*x2-27477/10000*x+5437/1000 多项式为Polyfit的第三个参数是多项式的阶数。,多项式积分：功能：求多项式积分
32、调用格式：py=poly_itg(p)p:被积多项式的系数 py：求积后多项式的系数 poly_itg.m function py=poly_itg(p)n=length(p);py=p.*n:-1:1.(-1),0不包括最后一项积分常数,多项式微分：Polyder:求多项式一阶导数的系数。调用格式为:b=polyder(c)c为多项式y的系数，b是微分后的系数，其值为：,两个多项式的和与差：命令poly_add:求两个多项式的和,其调用格式为:c=poly_add(a,b)多项式a减去b，可表示为:c=poly_add(a,-b),功能：两个多项式相加 调用格式：b=poly_add(p1,
33、p2)b:求和后的系数数组poly_add.mfunction p3=poly_add(p1,p2)n1=length(p1);n2=length(p2);if n1=n2 p3=p1+p2;endif n1n2 p3=p1+zeros(1,n1-n2),p2;endif n1n2 p3=zeros(1,n2-n1),p1+p2;end,m阶多项式与n阶多项式的乘积是d=m+n阶的多项式:计算 系数的MATLAB命令是:c=conv(a,b)多项式 除多项式 的除法满足:其中 是商,是除法的余数。多项式 和 可由命令deconv算出。例：q,r=deconv(a,b),例 a=2,-5,6,-
34、1,9;b=3,-90,-18;c=conv(a,b)c=6-195 432-453 9-792-162 q,r=deconv(c,b)q=2-5 6-1 9r=0 0 0 0 0 0 0 poly2sym(c)ans=6*x6-195*x5+432*x4-453*x3+9*x2-792*x-162,4.4 数据插值4.4.2 Lagrange插值,方法介绍 对给定的n个插值点 及对应的函数值，利用构造的n-1次Lagrange插值多项式，则对插值区间内任意x的函数值y可通过下式求的：MATLAB实现,function y=lagrange(x0,y0,x)ii=1:length(x0);y=
35、zeros(size(x);for i=ii ij=find(ii=i);y1=1;for j=1:length(ij),y1=y1.*(x-x0(ij(j);end y=y+y1*y0(i)/prod(x0(i)-x0(ij);end算例：给出f(x)=ln(x)的数值表，用Lagrange计算ln(0.54)的近似值。x=0.4:0.1:0.8；y=-0.916291,-0.693147,-0.510826,-0.356675,-0.223144;lagrange(x,y,0.54,0.55,0.78)ans=-0.6161-0.5978-0.2484（精确解-0.616143),4.4.
36、3 Hermite插值,方法介绍 不少实际问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式就是Hermite插值多项式。下面只讨论函数值与一阶导数值个数相等且已知的情况。已知n个插值点 及对应的函数值 和一阶导数值。则对插值区间内任意x的函数值y的Hermite插值公式：,MATLAB实现%hermite.mfunction y=hermite(x0,y0,y1,x)n=length(x0);m=length(x);for k=1:m yy=0.0;for i=1:n h=1.0;a=0.0;for j=1:n if j=i h=h*(x
37、(k)-x0(j)/(x0(i)-x0(j)2;a=1/(x0(i)-x0(j)+a;end end yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i);end y(k)=yy;end,算例:对给定数据,试构造Hermite多项式求出sin0.34的近似值。x0=0.3,0.32,0.35;y0=0.29552,0.31457,0.34290;y1=0.95534,0.94924,0.93937;y=hermite(x0,y0,y1,0.34)y=0.3335 sin(0.34)与精确值比较ans=0.3335,x=0.3:0.005:0.35;y=hermi
38、te(x0,y0,y1,x);plot(x,y)y2=sin(x);hold on plot(x,y2,-r),4.4.4 Runge现象,问题的提出：根据区间a,b上给出的节点做插值多项式p(x)的近似值，一般总认为p(x)的次数越高则逼近f(x)的精度就越好，但事实并非如此。反例：在区间-5,5上的各阶导数存在，但在此区间上取n个节点所构成的Lagrange插值多项式在全区间内并非都收敛。取n=10,用Lagrange插值法进行插值计算。,x=-5:1:5;y=1./(1+x.2);x0=-5:0.1:5;y0=lagrange(x,y,x0);y1=1./(1+x0.2);绘制图形 pl
39、ot(x0,y0,-r)插值曲线 hold on plot(x0,y1,-b)原曲线为解决Rung问题，引入分段插值。,算法分析：所谓分段插值就是通过插值点用折线或低次曲线连接起来逼近原曲线。MATLAB实现 可调用内部函数。命令1 interp1 功能：一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。格式1 yi=interp1(x,Y,xi)%返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。算例 对于t，beta、alpha分别有两组数
40、据与之对应，用分段线性插值法计算当t=321,440,571时beta、alpha的值。,4.4.5 分段插值,temp=300,400,500,600;beta=1000*3.33,2.50,2.00,1.67;alpha=10000*0.2128,0.3605,0.5324,0.7190;ti=321,400,571;propty=interp1(temp,beta,alpha,ti);propty=interp1(temp,beta,alpha,ti,linear);ti,proptyans=1.0e+003*0.3210 3.1557 2.4382 0.4000 2.5000 3.60
41、50 0.5710 1.7657 6.6489,格式2 yi=interp1(Y,xi)%假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。格式3 yi=interp1(x,Y,xi,method)%用指定的算法计算插值：nearest：最近邻点插值，直接完成计算；linear：线性插值（缺省方式），直接完成计算；spline：三次样条函数插值。cubic：分段三次Hermite插值。其它，如v5cubic。对于超出x范围的xi的分量，使用方法nearest、linear、v5cubic的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。yi=in
42、terp1(x,Y,xi,method,extrap)yi=interp1(x,Y,xi,method,extrapval)%确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。,算例 year=1900:10:2010;product=75.995,91.972,105.711,123.203,131.669,.150.697,179.323,203.212,226.505,249.633,256.344,267.893;p1995=interp1(year,product,1995)p1995=252.9885 x=1900:1:2010;y=interp1(ye
43、ar,product,x,cubic);plot(year,product,o,x,y),例:已知的数据点来自函数根据生成的数据进行插值处理，得出较平滑的曲线直接生成数据。x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);plot(x,y,x,y,o),调用interp1()函数：x1=0:.02:1;y0=(x1.2-3*x1+5).*exp(-5*x1).*sin(x1);y1=interp1(x,y,x1);y2=interp1(x,y,x1,cubic);y3=interp1(x,y,x1,spline);y4=interp1(x,y,x1,near
44、est);plot(x1,y1,y2,y3,y4,:,x,y,o,x1,y0)误差分析 max(abs(y0(1:49)-y2(1:49),max(abs(y0-y3),max(abs(y0-y4)ans=0.0177 0.0086 0.1598,x0=-1+2*0:10/10;y0=1./(1+25*x0.2);x=-1:.01:1;y=lagrange(x0,y0,x);%Lagrange 插值 ya=1./(1+25*x.2);plot(x,ya,x,y,:),例,y1=interp1(x0,y0,x,cubic);y2=interp1(x0,y0,x,spline);plot(x,ya
45、,x,y1,:,x,y2,-),命令2 interp2 功能 二维数据内插值（表格查找）格式1 ZI=interp2(X,Y,Z,XI,YI)%返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回NaN。格式2 ZI=interp2(Z,XI,YI)%缺省地，X=1:n、Y=1:m，其中m,n=size(Z)。再按第一种情形进行计算。格式3 ZI=interp2(X,Y,Z,XI,YI,method)%用指定的算法method计算二维插值：
46、linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。,算例：years=1950:10:1990;service=10:10:30;wage=150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243;w=interp2(service,years,wage,15,1975)w=190.6288,例 x,y=meshgrid(-3:.6:3,-2:
47、.4:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);surf(x,y,z),axis(-3,3,-2,2,-0.7,1.5),选较密的插值点，用默认的线性插值算法进行插值 x1,y1=meshgrid(-3:.2:3,-2:.2:2);z1=interp2(x,y,z,x1,y1);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5),立方和样条插值：z1=interp2(x,y,z,x1,y1,cubic);z2=interp2(x,y,z,x1,y1,spline);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.
48、5)figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5),算法误差的比较 z=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z-z1),axis(-3,3,-2,2,0,0.08)figure;surf(x1,y1,abs(z-z2),axis(-3,3,-2,2,0,0.025),二维一般分布数据的插值,功能:可对非网格数据进行插值格式:z=griddata(x0,y0,z0,x,y,method)v4:MATLAB4.0提供的插值算法，公认效果较好；linear：双线性插值算法（缺省算法）；ne
49、arest：最临近插值；spline：三次样条插值；cubic：双三次插值。例：在x为-3，3，y为2，2矩形区域随机选择一组坐标，用 v4 与cubic插值法进行处理，并对误差进行比较。,x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);x1,y1=meshgrid(-3:.2:3,-2:.2:2);z1=griddata(x,y,z,x1,y1,cubic);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)z2=griddata(x,y,z,x1,y1,v4);figur
50、e;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5),误差分析 z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.15)figure;surf(x1,y1,abs(z0-z2),axis(-3,3,-2,2,0,0.15),例：在x为-3，3，y为-2，2矩形区域随机选择一组坐标中，对分布不均匀数据，进行插值分析。x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);%生
