1、MATLAB语言及应用,PART3 数值计算,曾然 杭州电子科技大学 通信工程学院,数值计算内容:,数值计算包括:多项式运算, 拟合与插值，数值微积分，线性方程组求解, 数据分析, 信号处理和系统分析.,1 多项式运算,多项式运算是数学中最基本的运算之一.多项式一般可以表示为:f(x)=a0xn+a1xn-1+a2xn-2+an-1x+an 对于这种表示形式，很容易用一个行向量表示，即：T=a0,a1,a2,an-1,an 在MATLAB中, 多项式正是用这样一个行向量表示的,系数是按降序排列的,幂,多项式构造,多项式可以直接用向量表示, 因此, 构造多项式最简单的方法是直接输入向量.,【例】
2、 直接输入向量构造 f(x)=2x5+5x4+4x2+x+4,T=2, 5, 0, 4, 1, 4; fx=poly2sym(T),函数poly2sym是符号工具箱中的函数,在用此种方式构造多项式时, 必须把多项式各项的系数写完整, 而不管此项的系数是否为0,fx = 2*x5+5*x4+4*x2+x+4,r=1,2,3,4; T1=poly(r); y=poly2sym(T1) y_class=class(y) %显示出y为符号变量,【例】 用多项式的根构造多项式,根为r=1,2,3,4,T1 =1 -10 35 -50 24y =x4-10*x3+35*x2-50*x+24y_class
3、=sym,多项式的运算,多项式的运算主要包括多项式的四则运算,导数运算,估值运算,求根运算以及多项式的拟合、插值等,1 多项式的四则运算,多项式四则运算主要是多项式的加,减,乘,除.其中,多项式的加减运算要求两个相加,减的多项式的大小必须相等,此时加,减才有效. 当两个相加,减的多项式阶次不同时必须通过补0使其相同.,加/减-+/-,乘/除-conv, deconv,T, r=deconv(T1,T3),商多项式,余式,T3为分母,T1=2, 5, 0, 4, 1, 4;T2=0, 0, 5, 1, 3, 2;T3=5, 1, 3, 2; % 重写T2，因除法运算中分母多项式第一个系数不能为0
4、（为了与T1做加法运算的T2前两个系数补了0，不能再作除法运算）T=T1+T2; % 必须是同维的才能相加T_add=poly2sym(T)T=T1-T2;T_sub=poly2sym(T) % 注意：加减运算只要对系数操作即可，但乘除运算需要用conv和deconv指令T=conv(T1,T2); % 乘法不要求同维T_mul=poly2sym(T)A_coe, A_r=deconv(T1,T3);T_coe=poly2sym(A_coe)T_rem=poly2sym(A_r),【例】 多项式的加减乘除运算,f1(x)=2x5+5x4+4x2+x+4, f2(x)=5x3+x2+3x+2,【
5、例】 多项式求值,求上式f1(x)在x=0.5处的函数值,T1=2,5,0,4,1,4;x=0.5;y=polyval(T1,x),y =5.8750,Polyval(p,X)求多项式的值指令： 当X为输入数组时，则对应求出每个元素作为多项式自变量求出的各个函数值作为输出数组,3 多项式的导数运算-polyder,【例】 求多项式f1(x)=2x5+5x4+4x2+x+4的导数,T1=2,5,0,4,1,4;h=polyder(T1);poly2sym(h),2 多项式求根-roots,【例】 求多项式f1(x)=2x5+5x4+4x2+x+4的根,T1=2, 5, 0, 4, 1, 4;ro
6、ot=roots(T1),root =-2.7709 0.5611 + 0.7840i0.5611 - 0.7840i-0.4257 + 0.7716i-0.4257 - 0.7716i,10*x4+20*x3+8*x+1,4 拟合和插值-polyfit, interp1,p = polyfit(x0, y0, n); p, s = polyfit(x0, y0, n);,x0,y0-给定数据对 n-拟合出的多项式次数 p-多项式向量 s-偏差信息,yi=interp1(x0, y0, xi, cubic);,xi, yi-得到的新的插值点 cubic-插值方式,插值方式:,linear-线性
7、插值 nearest-最近插值 cubic-三次多项式插值 spline-样条插值,拟合和插值的区别:曲线拟合是研究如何寻找“平滑”曲线最好地表现带噪声的“测量数据”,但并不要求拟合曲线通过这些“测试数据”点。插值是在认定所给数据完全正确的情况下,研究如何平滑地估算出“基准数据”之间其他点的函数值,因此插值一定通过“基准数据”。,【拟合的例子】 由下列数据对(x0,y0)求三次和六次拟合多项式并绘图.x0=0:0.1:1;y0=-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22;,通俗来讲, 拟合就是由已知点得到一条曲线, 使该曲
8、线最能反映点所代表的规律.比如做欧姆定理的实验的时候,由于实验中存在误差,最后拟合得到的曲线是一条直线,而且肯定只有部分点落在拟合的直线上,但此时该直线和测试点的方差最小.由拟合直线的斜率就可以知道电阻的阻值.拟合是探测事物变化规律的办法. 插值就是根据函数上某些已知点(或实验数据),按一定规律(插值方法)寻求未知的点,比如已知一个常用对数y=log(x)表,是按照x=0.1:0.1:10制表的,如果按已知数据求y=log(2.897)就可以用插值得到.表制得越密,插值越准确,x0=0:0.1:1; y0=-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3
9、.45,5.35,9.22; n=3; % 设定拟合次数为3 p,s=polyfit(x0,y0,n); % 得到拟合多项式向量和相关偏差信息 xx=0:0.01:1; yy=polyval(p,xx); % 按拟合曲线计算采样值，以便画图 n1=6; % 设定拟合次数为6 p1,s1=polyfit(x0,y0,n1); yy1=polyval(p1,xx); plot(xx,yy,-b,xx,yy1,-m,x0,y0,.r,MarkerSize,20); y3=polyval(p,0.5); y6=polyval(p1,0.5); text(0.5, y3-0.3, n=3); text(
10、0.5, y6+0.2, n=6);,注意：拟合只能在给定数据所限定的区间内使用, 不要任意向外拓展,y=51.48x3 -77.74x2+35.06x-0.20,y=348.21x6-1060.23x5+1297.68x4-758.90x3+ 181.00x2 + 1.00x+0.48,【例】 按上例所给数据研究插值, 并观察插值和拟合的区别.,x0=0:0.1:1; y0=-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22; xi=0:0.02:1; yi=interp1( x0, y0, xi, cubic ); subpl
11、ot(1,2,1); plot(xi, yi, -b, x0, y0, .r, MarkerSize, 20); xlabel(x), ylabel(y); p,s=polyfit(x0,y0,3); xx=0:0.01:1; yy=polyval(p,xx); subplot(1,2,2); plot(xx, yy, -r, x0, y0, .r, MarkerSize, 20); xlabel(x); ylabel(y);,插值,拟合,2 数据分析函数,(1) 若输入宗量x是向量,那么不管是行向量还是列向量,运算是对整个向量进行的 (2) 若输入宗量x是二维数组，那么指令运算是按列进行的
12、见下列例子,MATLAB在进行数据分析时的约定:,(1) 随机数发生器和统计分析指令,【例】 基本统计示例,randn(state,0); A=randn(1000,4); AMAX=max(A), AMIN=min(A) %在+3,-3附近 AMED=median(A), AMEAN=mean(A) %接近0 ASTD=std(A) %接近1,已知A=1,2,3;4,5,6;7,8,9, 演示max, min, median, mean, std, var,各列最大,最小值：max(A), min(A) 各行最大,最小值：max(A), min(A) 各列中间值：median(A) 各列平均
13、值：mean(A) 各列方差,标准差： var(A), std(A),A= 1, 2, 3;4, 5, 6;7, 8, 9 ,7 8 9,(2) 差分和累计指令,(2) 差分和累计指令,注： diff:当X是向量时，dx=X(2:n)-X(1:n-1)；当X是矩阵时，dx=X(2:n,:)-(1:n-1, :)即按列求差分。注意diff(A)的长度比X的长度短少一个元素。例：A=1 4 2; 3 6 9; 3 2 5, diff(A), diff(A,1,1), diff(A,1,2) gradient:输入z为数组时， DZx,DZy给出与z同样大小的数组， DZx的每行给出z相应行元素之间
14、的“梯度”，即如第一行DZx(1,1)=Z(1,2)-Z(1,1), DZx(1,end)=Z(1,end)-Z(1,end-1), DZx(1,2:end)=1/2*Z(1,3:end)-Z(1,1:end-2)； DZy的每列给出z相应列元素之间的“梯度”。例：Z=1 2 5 4;3 2 4 5;3 10 5 6;4 8 6 3, DZx,DZy=gradient(Z) cumsum与cumtrapz的区别（矩形法、梯形法） trapz(x,y)给出采样点(x,y)所连接折线下的积分面积（为单值）；S=cumtrapz(x,y)为累计积分，即输出为一个与y大小相同的数组，每个元素S(k)为
15、对应该位置的自变量区间x(1)x(k)上的函数y积分值,diff, trapz, cumtrapz指令单参数输入的演示,a=1,2,3,4,5; b=a, a+1, a+2,a+3,diff(a)=diff(a,1)=1,1,1,1 diff(b)=1,1,1,1 ; 1,1,1,1 ; 1,1,1,1 ; 1,1,1,1 trapz(b)=12 16 20 24 cumtrapz(b),0 0 0 0 1.50 2.50 3.50 4.50 4.00 6.00 8.00 10.00 7.50 10.50 13.50 16.50 12.00 16.00 20.00 24.00,1 2 3 4
16、2 3 4 5 3 4 5 6 4 5 6 7 5 6 7 8,【例】(非累计的指令sum和prod:给出单值而非数组) 求1+2+3+100以及50!,x=1:1:100; sum_x=sum(x);,sum_x =5050,a=1:1:50; prod_a=prod(a),prod_a =3.0414e+064,cumsum, cumprod的用法,分别是求向量的累计和与累计乘积.如果a=1:1:n, 则 cumsum(a)=1, 1+2, 1+2+3, ., 1+2+3+n cumprod(a)=1, 1*2, 1*2*3, , 1*2*3*n,【例】 求f(x)=3x2在区间0,2的积
17、分,dt=0.001; t=0:dt:2; y=3*t.2; s1=dt*sum(y) s2=dt*trapz(y) s=dt*cumsum(y); s3=s(end) s=dt*cumtrapz(y); s4=s(end),matlab还有更精良的积分指令: quad, quadl,s1 =8.0060,s2 =8.000001,s3 =8.0060,s4 =8.000001,3 数值微积分,（1）数值积分,数值积分有闭型和开型算法之分两者的区别在于：是否需要计算积分区间端点处的函数值,matlab中仅提供闭型数值积分指令:,q=quad(fun, a, b, tol)-采用递推自适应Sim
18、pson法计算积分（低阶法数值积分）,q=quadl(fun, a, b, tol)-采用递推自适应Lobatto法求数值积分（高阶法数值积分）,SS=dblquad(fun, inmin, inmax, outmin, outmax,tol) -二重积分指令,fun-被积函数, a, b -积分上下限 tol-绝对误差控制量，缺省时积分的绝对精度为1x10-6,例4.4.3-1求f(x)=3exp(x2)在区间0,2的积分,Format long y=3*exp(x.2); a=0; b=2; q1=quad(y,a,b) q2=quadl(y,a,b) q3=quadl(y,a,b,1e-
19、8),q1 =49.35788332040026 q2 =49.35788329741851 q3 =49.35788329652209,（2） 常微分方程的数值解,matlab为解常微分方程初值问题提供了一组配套齐全,结构严整的指令, 包括: ode45, ode23, ode113, ode23t, ode15s, ode23s, ode23tb。最常用的是ode45，介绍其基本使用方法.,ode45使用方法:,x, y=ode45(file_name, xspan, y0),输出量：x-微分区间的点系列,y-原函数在微分区间点系列上的函数值,输入量：file_name-文档名,保存的M文
20、件名（该函数文件的输出必须是待解函数的一阶导数，即y=f(y,x)）,xspan-微分区间（一般为二元向量x0,xf）,y0-初始值,【例】 解微分方程：y=3x2, y(0)=1,微分区间0,5,% M文件,保存为dy_file.m function dy=f(x,y) dy=3*x2;,y0=1; xspan=0,5; x, y=ode45(dy_file,xspan,y0); plot(x,y,r), xlabel(x), ylabel(y),（3）求函数极值点,许多科学研究和工程计算问题都可以归结为一个极值问题,如能量最小，时间最短,最佳拟合,最短路径等.又因为f(x)的极小值问题等价
21、于-f(x)的极大值问题,所以matlab只有处理极小值的指令,确切地说,这里讨论的只是“局域极值”问题.“全域最小”问题要复杂得多.至今没有一个系统性的方法可求解一般的“全域最小”问题.对于一元二元函数,作图观察对确定全域最小值有很好的应用价值.但更多元的函数,就很难使用作图法.,matlab求函数极值的三条指令,x, fval, exitflag, output=fminbnd(fun,x1,x2,opt)-求一元函数在区间(x1,x2)中极小值的指令,x, fval, exitflag, output=fminsearch(fun,x0,opt)-单纯形法求多元函数极值点的指令*,x,f
22、val,exitflag,output,grad,hessian=fminunc(fun,x0,opt)-拟牛顿法求多元函数极值点的指令*,【例】 求函数y=2x-10x的最小值,(1) 函数采用字符串表达式 func=2x-10*x; x1=0;x2=8; % 通过观察, 发现最小值在此区间 Xmin,Ymin=fminbnd(func,x1,x2),Xmin =3.8507 Ymin =-24.0800,(2) 函数采用M文件形式 % M文件内容，保存为func.m function y=f(x) y=2x-10*x; % M文件内容 x1=0; x2=8;Xmin,Ymin=fminbn
23、d(func,x1,x2),xx=0:0.01:8; yy=2.xx-10.*xx; plot(xx,yy),4 信号处理和系统分析举例数字滤波,y=filter(B, A, x)-运用数字滤波对输入信号x滤波,滤波器的作用:通过对采样数据信号进行数学运算处理来达到频域滤波的目的.,例 设计一个低通滤波器,从受噪声干扰的多频率混合信号x(t) 中获取10HZ的信号. x(t)=sin(2*10*t)+cos(2*100*t) +n(t).,randn(state,1); ws=1000; % 采样频率 t=0:1/ws:0.4; % 产生带噪声的多频率信号 x=sin(2*pi*10*t)+cos(2*pi*100*t)+0.2*randn(size(t); wn=ws/2;fc=30; % Nyquest频率wn,绝对截止频率fc B,A=butter(10,30/wn); % 产生截止频率为fc/wn的10阶Butterworth低通滤波器分子 、分母多项式B和A，作为下面滤波处理的输入 y=filter(B, A, x); % 滤波处理 plot(t,x,b-,t,y,r.,MarkerSize,10);,
