《《数值计算》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《数值计算》PPT课件.ppt（50页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第四章第四章 数值计算数值计算1数值计算内容数值计算内容:数值计算包括数值计算包括:多项式运算多项式运算,线性方程组求解线性方程组求解,矩阵特征值问题矩阵特征值问题的解的解,卷积卷积,数据分析数据分析,泛函指令的使用泛函指令的使用,信号处理和系统分析信号处理和系统分析.4.1 多项式运算多项式运算多项式运算是数学中最基本的运算之一多项式运算是数学中最基本的运算之一.多项式一般可多项式一般可以表示为以表示为:f(x)=a0 xn+a1xn-1+a2xn-2+an-1x+an对于这种表示形式，很容易用一个行向量表示，即：对于这种表示形式，很容易用一个行向量表示，即：T=a0,a1,a2,an-1,
2、an在在MATLAB中中,多项式正是用这样一个行向量表示的多项式正是用这样一个行向量表示的,系数是按降序排列的系数是按降序排列的幂幂24.1.1 多项式构造多项式构造多项式可以直接用向量表示多项式可以直接用向量表示,因此因此,构造多项式最简单的方法是构造多项式最简单的方法是直接输入向量直接输入向量.例例4.1.1-1 直接输入向量构造直接输入向量构造 f(x)=2x5+5x4+4x2+x+4T=2,5,0,4,1,4;fx=poly2sym(T)函数函数poly2sym是符号工具箱中的函数是符号工具箱中的函数,在用此种方式在用此种方式构造多项式时构造多项式时,必须把多项式各项的系数写完整必须把
3、多项式各项的系数写完整,而不而不管此项的系数是否为管此项的系数是否为0fx=2*x5+5*x4+4*x2+x+43r=1,2,3,4;T1=poly(r);y=poly2sym(T1)y_class=class(y)例例4.1.1-2 用多项式的根构造多项式用多项式的根构造多项式,根为根为r=1,2,3,4T1=1 -10 35 -50 24y=x4-10*x3+35*x2-50*x+24y_class=sym44.1.2 多项式的运算多项式的运算多项式的运算主要包括多项式的四则运算多项式的运算主要包括多项式的四则运算,导数运算导数运算,估值运算估值运算,求根运算以及多项式的拟合等求根运算以及
4、多项式的拟合等1 多项式的四则运算多项式的四则运算多项式四则运算主要是多项式的加多项式四则运算主要是多项式的加,减减,乘乘,除除.其中其中,多多项式的加减运算要求两个相加项式的加减运算要求两个相加,减的多项式的大小必须减的多项式的大小必须相等相等,此时加此时加,减才有效减才有效.当两个相加当两个相加,减的多项式阶次减的多项式阶次不同时必须通过补不同时必须通过补0使其相同使其相同.加加/减减-+/-乘/除除-conv,deconvT=deconv(T1,T3)T,r=deconv(T1,T3)商多项式商多项式余式余式T3为分母为分母5T1=2,5,0,4,1,4;T2=0,0,5,1,3,2;T
5、3=5,1,3,2;%除法运算中分母多项式第一个系数不能为除法运算中分母多项式第一个系数不能为0T=T1+T2;%必须是同维的才能相加必须是同维的才能相加T_add=poly2sym(T)T=T1-T2;T_sub=poly2sym(T)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)例例4.1.1-3 多项式的加减乘除运算多项式的加减乘除运算f1(x)=2x5+5x4+4x2+x+4,f2(x)=5x3+x2+3x+26例
6、例4.1.1-4 多项式求值多项式求值,求上式求上式f1(x)在处的函数值在处的函数值T1=2,5,0,4,1,4;x=0.5;y=polyval(T1,x)y=73 多项式的导数运算多项式的导数运算-polyder例例4.1.1-6 求多项式求多项式f1(x)=2x5+5x4+4x2+x+4的导数的导数 T1=2,5,0,4,1,4;h=polyder(T1);poly2sym(h)2 多项式求根多项式求根-roots例例4.1.1-5 求多项式求多项式f1(x)=2x5+5x4+4x2+x+4的根的根 T1=2,5,0,4,1,4;root=roots(T1);root=-2.7709 1
7、0*x4+20*x3+8*x+184 拟合和插值拟合和插值-polyfit,interp1例例4.1.1-7 对下列数据对对下列数据对(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;p=polyfit(x0,y0,n);p,s=polyfit(x0,y0,n);x0,y0-给定数据对给定数据对n-拟合出的多项式次数拟合出的多项式次数p-多项式向量多项式向量s-偏差信息偏差信息yi=interp1(x0,y0,xi,cubic);xi,yi-
8、得到的新的插值点得到的新的插值点cubic-插值方式插值方式9x0=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;n=3;%设定拟合次数为设定拟合次数为3p,s=polyfit(x0,y0,n);%得到拟合多项式向量和相关偏差信得到拟合多项式向量和相关偏差信息息 xx=0:0.01:1;yy=polyval(p,xx);%按拟合曲线计算采样值按拟合曲线计算采样值n1=6;%设定拟合次数为设定拟合次数为6p1,s1=polyfit(x0,y0,n1);yy1=polyval(p1,xx);plot(xx,
9、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(0.5,y6+0.2,n=6);10326543+211拟合只能在给定数据所限定的区间内使用拟合只能在给定数据所限定的区间内使用,不要任意向外拓展不要任意向外拓展例例4.1.1-8 按上例所给数据研究插值按上例所给数据研究插值,并观察插值和拟合的区别并观察插值和拟合的区别.x0=0:0.1:1;y0=-0.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.4
10、5,5.35,9.22;xi=0:0.02:1;yi=interp1(x0,y0,xi,cubic);subplot(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);12插值插值拟合拟合13插值方式插值方式:linear-线性插值线性插值nearest-最近插值最近插值c
11、ubic-三次多项式插值三次多项式插值spline-样条插值样条插值pchip-分段三次分段三次Hermite插值插值版的三次多项式插值版的三次多项式插值.拟合和插值的区别拟合和插值的区别:曲线拟合是研究如何寻找曲线拟合是研究如何寻找“平滑平滑”曲线最好地表现带噪声的曲线最好地表现带噪声的“测量数据测量数据”,但并不要求拟但并不要求拟合曲线通过这些合曲线通过这些“测试数据测试数据”点点.插值是在认定所给数插值是在认定所给数据完全正确的情况下据完全正确的情况下,研究如何平滑地估算出研究如何平滑地估算出“基准数基准数据据”之间其他点的函数值之间其他点的函数值,因此插值一定通过因此插值一定通过“基准
12、数基准数据据”14通俗地讲通俗地讲,拟合就是由已知点得到一条曲线拟合就是由已知点得到一条曲线,使该曲线使该曲线最能反映点所代表的规律最能反映点所代表的规律.比如做欧姆定理的实验的时比如做欧姆定理的实验的时候候,由于实验中存在误差由于实验中存在误差,最后拟合得到的曲线是一条最后拟合得到的曲线是一条直线直线,而且肯定只有部分点落在拟合的直线上而且肯定只有部分点落在拟合的直线上,但此时但此时该直线和测试点的方差最小该直线和测试点的方差最小.由拟合直线的斜率就可以由拟合直线的斜率就可以知道电阻的阻值知道电阻的阻值.拟合是探测事物变化规律的办法拟合是探测事物变化规律的办法.插值就是根据函数上某些已知点插
13、值就是根据函数上某些已知点(或实验数据或实验数据),按一定按一定规律规律(插值方法插值方法)寻求未知的点寻求未知的点,比如已知一个常用对数比如已知一个常用对数y=log(x)表表,是按照是按照x=0.1:0.1:10制表的制表的,如果按已知数如果按已知数据求据求y=log(2.897)就可以用插值得到就可以用插值得到.表制得越密表制得越密,插插值越准确值越准确15例例4.1.1-9 已知通过实验得到电压和电流的数值如下：已知通过实验得到电压和电流的数值如下：i0=0:0.5:10;v0=0,98,202,306,395,506,592,708,789,890,1009,1108,1186,13
14、10,1391,1510,1601,1716,1782,1920,2014;按一次按一次,二次二次,三次多项式对数值关系进行拟合三次多项式对数值关系进行拟合说明：说明：为保证较好的拟合效果，多项式阶数要取得适当，过为保证较好的拟合效果，多项式阶数要取得适当，过低，残差较大，过高，拟合模型将包含噪声影响，通低，残差较大，过高，拟合模型将包含噪声影响，通常要保证拟合阶数小于数据对的数目。常要保证拟合阶数小于数据对的数目。16i0=0:0.5:10;v0=0,98,202,306,395,506,592,708,789,890,1009,1108,.1186,1310,1391,1510,1612,
15、1716,1782,1920,2014;n1=1;p1,s=polyfit(i0,v0,n1);n2=2;p2,s=polyfit(i0,v0,n2);n3=3;p3,s=polyfit(i0,v0,n3);ii=0:0.1:10;v1=polyval(p1,ii);v2=polyval(p2,ii);v3=polyval(p3,ii);subplot(1,3,1),plot(ii,v1,-b,i0,v0,.r,MarkerSize,8),xlabel(一次拟合一次拟合);subplot(1,3,2),plot(ii,v2,-b,i0,v0,.r,MarkerSize,8),xlabel(二次
16、拟合二次拟合);subplot(1,3,3),plot(ii,v3,-b,i0,v0,.r,MarkerSize,8),xlabel(三次拟合三次拟合);p1=201.00,-2.89P2=0.31,197.90,2.02P3=0.03,-0.10,199.49,0.8517p1=201.00,-2.89P2=0.31,197.90,2.02P3=0.03,-0.10,199.49,0.8518对于含对于含n个未知数的个未知数的n个方程构成的方程组个方程构成的方程组Ax=b,在线性代数教在线性代数教科书中科书中,最常介绍的解法有最常介绍的解法有:Cramer法法;逆阵法逆阵法,即即x=A-1b
17、;高斯消高斯消元法元法;LU法法 对于维数不高对于维数不高,条件数不大的矩阵条件数不大的矩阵,以上四种所得结果不会有明以上四种所得结果不会有明显的差别显的差别.但前三种解法的更多的意义在理论上但前三种解法的更多的意义在理论上,而不在实际的而不在实际的数值计算上数值计算上.在在MATLAB中中,方程采用方程采用LU法求解法求解,并且出于算法并且出于算法稳定性的考虑稳定性的考虑,行列式和逆的计算也都是在行列式和逆的计算也都是在LU分解的基础上进分解的基础上进行的行的MATLAB中的四条指令中的四条指令:L,U,P=lu(A)-矩阵的矩阵的LU分解分解,使使LU=PA.多种格式多种格式.det(A)
18、-求矩阵求矩阵A的行列式的行列式,它由它由detA=detUuii 算得算得inv(A)-求矩阵求矩阵A的逆的逆,由由A-1=U-1L-1P算得算得x=Ab-除法指令求方程的解除法指令求方程的解(对恰定方程对恰定方程,采用采用LU法执行法执行)4.2.1 方程组求解方程组求解4.2 线性方程组的解线性方程组的解19对于方程组对于方程组Ax=b,采用采用x=Ab计算，如果方程组为计算，如果方程组为yC=d,要使用右除，即指令为要使用右除，即指令为y=d/CAx=bxA=byC=dx=Abx=b/Ay=d/C例例4.2.1-1 解下列方程组解下列方程组2x1+2x2+3x3=34x1+7x2+7x
19、3=1-2x1+4x2+5x3=-7A=2,2,3;4,7,7;-2,4,5,b=3;1;-7x=Ab,y=b/A求解方程时求解方程时,尽量不要使用指令尽量不要使用指令inv(A)*b进行进行,它不仅计算速它不仅计算速度没有除法快度没有除法快,而且计算精度也没除法高而且计算精度也没除法高.此外此外,除法还适用于除法还适用于“超定超定”和和“欠定欠定”方程方程x=y=20方程组的解的三种情况方程组的解的三种情况:对于方程对于方程Ax=b,A为为Amn矩阵矩阵,有三种情况有三种情况:当当m=n时时,此方程成为此方程成为恰定恰定方程方程,求解精确解求解精确解 当当mn时时,此方程成为此方程成为“超定
20、超定”方程方程,寻求寻求最小二乘解最小二乘解(直线拟直线拟合合)当当mn时此时无解时此时无解,但实用中可以求最小二乘解即但实用中可以求最小二乘解即:方程解方程解(AA)x=A b x=(AA)-1Ab-求逆法求逆法 x=Ab-matlab用最小二乘法找一个近似解用最小二乘法找一个近似解.x1+2x2=1 2x1+3x2=23x1+4x2=3例例4.2.1-3 求下列超定方程的解求下列超定方程的解 求逆法：求逆法：x=inv(A*A)*A*b 使用使用Ab 的方式的方式A=1,2;2,3;3,4;b=1;2;3;x1=inv(A*A)*A*b;x2=Abx1=x2=022 3)3)欠定方程欠定方
21、程当方程数少于未知量个数时当方程数少于未知量个数时,即不定情况即不定情况,有无穷多个解存在有无穷多个解存在.matlab可求出两个解：可求出两个解：a.用除法求的解用除法求的解x是具有最多零元素的解是具有最多零元素的解b.具有最小长度或范数的解具有最小长度或范数的解,这个解是基于伪逆这个解是基于伪逆pinv求得的求得的.例求下列欠定方程组的解例求下列欠定方程组的解x1+2x2+3x3=12x1+3x2+4x3=2A=1,2,3;2,3,4;b=1;2;x1=Ab;x2=pinv(A)*bx1=1 0 0 x2=234.2.2 奇异值分解和矩阵结构奇异值分解和矩阵结构相应的相应的MATLAB指令
22、如下指令如下:U,S,V=svd(A)-矩阵的奇异值分解三对组阵，使矩阵的奇异值分解三对组阵，使A=USVTnorm(A,flag)-计算矩阵计算矩阵A的范数，范数类型由的范数，范数类型由flag指定指定.cond(A)-计算矩阵计算矩阵A的条件数的条件数.pinv(A)-求矩阵的广义逆求矩阵的广义逆flag-1,2,inf,fro244.2.3 线性二乘问题的解线性二乘问题的解对于超定方程对于超定方程,求其最小二乘解有三种方法：求其最小二乘解有三种方法：正则方程法得解正则方程法得解x=(A A)-1A b,数值精度差数值精度差广义逆法得解广义逆法得解x=A+b,所得解可靠所得解可靠.用矩阵除
23、法得解用矩阵除法得解x=Ab.可靠性稍差，但速度快可靠性稍差，但速度快.例例4.2.3-1 对于超定方程对于超定方程y=Ax,进行三种解法比较进行三种解法比较.(1)生成矩阵生成矩阵A及及y,并用三种方法求解并用三种方法求解A=gallery(5),A(:,1)=;y=1.7,7.5,6.3,0.83,-0.082;x=inv(A*A)*A*y,xx=pinv(A)*y,xxx=Ayx=xx=xxx=0由于由于阵缺一个列秩阵缺一个列秩,除法给出的最小二乘除法给出的最小二乘解就只有三个非零元解就只有三个非零元素素.它是只有三个非零它是只有三个非零元素的所有最小二乘元素的所有最小二乘解中范数最小的
24、解中范数最小的Gallery-MATLAB 设置设置的特殊矩阵的特殊矩阵.25(2)计算三个解的范数计算三个解的范数nx=norm(x),nxx=norm(xx),nxxx=norm(xxx)nx=nxx=nxxx=用广义逆所得的最小二乘范数最小(3)比较三种解法的方程误差比较三种解法的方程误差e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx)e=ee=eee=误差的平方根最大，精度最差误差的平方根最大，精度最差264.2.4 特征值分解和矩阵函数特征值分解和矩阵函数涉及矩阵特征值分解的常用的涉及矩阵特征值分解的常用的MATLAB指令如下指令如下:d
25、=eig(A)-计算矩阵计算矩阵A的特征值的特征值,以向量形式存放以向量形式存放V,D=eig(A)-计算矩阵的特征向量计算矩阵的特征向量V和特征值对角阵和特征值对角阵D,使使A*V=V*DVR,DR=cdf2rdf(VC,DC)-把复数对角形转化成实数块对角形把复数对角形转化成实数块对角形.VC,DC=rsf2csf(VR,DR)-把实数块对角形转换成复数对角形把实数块对角形转换成复数对角形V,J=jordan(A)-jordan分解分解,使使A*V=V*J,J是对角阵是对角阵c=condeig(A)-计算各特征值的条件数计算各特征值的条件数V,D,C=condeig(A)-相当于相当于V,
26、D=eig(A)和和c=condeig(A)两条指令两条指令27例例4.2.4-1 验证验证det(A)=iA=2,-1,-1;3,4,-2;3,-2,4det(A)d=eig(A)ans=60d=6.0000 d(1)*d(2)*d(3)ans=行列式行列式特征值特征值%矩阵行列式等于其所有特征值的积矩阵行列式等于其所有特征值的积284.3 数据分析函数数据分析函数(1)若输入宗量若输入宗量x是向量是向量,那么不管是行向量还是列向量那么不管是行向量还是列向量,运算是运算是对整个向量进行的对整个向量进行的(2)若输入宗量若输入宗量x是二维数组，那么指令运算是按列进行的是二维数组，那么指令运算是
27、按列进行的MATLAB在进行数据分析时的约定在进行数据分析时的约定:294.3.1 随机数发生器和统计分析指令随机数发生器和统计分析指令函数函数含义含义rand(m,n)产生产生mxn维的维的0,1区间均匀分布随机数组区间均匀分布随机数组randn(m,n)产生产生mxn维的均值为维的均值为0标准差为标准差为1的正态分布随的正态分布随机数组机数组min(x)对矩阵各列分别求最小值对矩阵各列分别求最小值max(x)对矩阵各列分别求最大值对矩阵各列分别求最大值median(x)对矩阵各列分别求中位数对矩阵各列分别求中位数mean(x)对矩阵各列分别求均值对矩阵各列分别求均值std(x)对矩阵各列分
28、别求标准差对矩阵各列分别求标准差var(x)对矩阵各列分别求方差对矩阵各列分别求方差C=conv(x)给出矩阵各列的协方差阵给出矩阵各列的协方差阵P=corrcoef(x)给出矩阵各列间的相关系数给出矩阵各列间的相关系数30例例4.3-1 基本统计示例基本统计示例randn(state,0);A=randn(1000,4);AMAX=max(A),AMIN=min(A)AMED=median(A);AMEAN=mean(A);ASTD=std(A)已知已知A=1,2,3;4,5,6;7,8,9,演示演示max,min,median,mean,std,var,conv,corrcoefvar,c
29、onv,corrcoef各列最大各列最大,最小值：最小值：max(A),min(A)各行最大各行最大,最小值：最小值：max(A),min(A)各列中间值：各列中间值：median(A)各列平均值：各列平均值：mean(A)各列方差各列方差,标准差：标准差：var(A),std(A)A=1,2,3;4,5,6;7,8,9 7 8 9314.3.2 差分和累计指令差分和累计指令 函数函数含义含义del2(U,hx,hy)五点五点Laplaciandiff(X,m,n)沿第沿第n维求维求m次差分次差分(二维：二维：n=1;列列,n=2,行行)DZx,DZy=gradient(z,dx,dy)对对Z
30、求求x,y方向梯度方向梯度prod(X,n)沿第沿第n维求乘积维求乘积sum(X,n)沿第沿第n维求和维求和trapz(x,Y,n)用梯形法沿第用梯形法沿第n维求函数关于自变量维求函数关于自变量x的积分的积分cumprod(X,n)沿第沿第n维求累计乘积维求累计乘积cumsum(X,n)沿第沿第n维求累计和维求累计和cumtrapz(x,Y,n)用梯形法沿第用梯形法沿第n维求函数维求函数Y关于关于x的累计积分的累计积分XS,KK=sort(X,n)沿第沿第n维对维对X元素按模增大排列元素按模增大排列32diff,trapz,cumtrapz指令的演示指令的演示a=1,2,3,4,5;b=a,a
31、+1,a+2,a+3diff(a)=diff(a,1)=1,1,1,1diff(b)=1,1,1,1;1,1,1,1;1,1,1,1;1,1,1,1trapz(b)=12 16 20 24cumtrapz(b)0 0 0 01 2 3 42 3 4 53 4 5 64 5 6 75 6 7 833例求例求1+2+3+100以及以及50!x=1:1:100;sum_x=sum(x);sum_x=5050a=1:1:50;prod_a=prod(a)prod_a=3.0414e+064cumsum,cumprod的用法的用法上两个函数分别是求向量的累计和和累计乘积上两个函数分别是求向量的累计和和累
32、计乘积.如果如果a=1:1:n,则则cumsum(a)=1,1+2,1+2+3,.,1+2+3+ncumprod(a)=1,1*2,1*2*3,1*2*3*n34例例4.3.2-1 求求f(x)=3x2在区间在区间0,2的积分的积分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,quda8s1=s2=s3=s4=354.4 MATLAB泛函指令泛函指令在在MATLAB中中,
33、凡以函数为输入宗量的指令凡以函数为输入宗量的指令,都被统称为泛函指令都被统称为泛函指令.最常见的最常见的有有:z=fzero(fun,x0)-求一元函数零点指令的最简单格式求一元函数零点指令的最简单格式x=fsolve(fun,x0)-解非线性方程组的最简单格式解非线性方程组的最简单格式x=fminbnd(fun,x1,x2)-求函数在区间求函数在区间(x1,x2)中极小值的指令中极小值的指令 最简格式最简格式x=fminsearch(fun,x0)-单纯形法求多元函数极值点指令的最单纯形法求多元函数极值点指令的最简格式简格式.x=fminunc(fun,x0)-拟牛顿法求多元函数极值点指令的
34、最简拟牛顿法求多元函数极值点指令的最简格式格式a=lsqnonlin(fun,a0)-解非线性最小二乘问题指令的最简格式解非线性最小二乘问题指令的最简格式.36q=quad(fun,a,b)-采用递推自适应采用递推自适应simpson法计算积分法计算积分q=quadl(fun,a,b)-采用递推自适应采用递推自适应Lobatto法求数值积分法求数值积分SS=dblquad(fun,inmin,inmax,outmin,outmax)-二重闭环数二重闭环数值积分值积分t,YY=ode45(fun,tspan,Y0)-采用采用4,5阶阶Runge-Kutta方程方程解算解算ODE初值问题初值问题指
35、令中被处理的函数指令中被处理的函数fun,可以取可以取:字符串表达式字符串表达式,内联函数内联函数,M函数函数文件的函数句柄文件的函数句柄374.4.1 求函数零点求函数零点对于任意函数对于任意函数f(x)=0,零点情况复杂零点情况复杂,没有一个通用解法没有一个通用解法.一般来说一般来说,零点的数值计算过程是零点的数值计算过程是:先猜测一个初始零点先猜测一个初始零点或该零点所在的区间或该零点所在的区间,然后通过计算然后通过计算,使猜测值不断精确使猜测值不断精确化化,或使猜测区间不断收缩或使猜测区间不断收缩,直到达到预定的精度直到达到预定的精度,终止计终止计算算.求函数的零点的步骤为求函数的零点
36、的步骤为:(1)利用利用matlab作图指令获取初步近似值作图指令获取初步近似值具体做法具体做法:先确定一个零点可能存在的自变量区间先确定一个零点可能存在的自变量区间,然后利用然后利用plot指令画出指令画出f(x)在该区间中的图形在该区间中的图形,观察观察f(x)与横轴的交点坐标与横轴的交点坐标,也可也可以用以用zoom对交点处局部放大再读数对交点处局部放大再读数.借助借助ginput指令获得更精确指令获得更精确的交点坐标值的交点坐标值38(2)利用利用fzero指令求精确解指令求精确解z=fzero(fun,x0)-求一元函数零点指令的最简格式求一元函数零点指令的最简格式z,f_z,exi
37、tflag=fzero(fun,x0,options,p1,p2,)-最完整格式最完整格式 说明：说明：由于由于fzero是根据函数是否越横轴来决定零点的是根据函数是否越横轴来决定零点的,因此本程序因此本程序无法确定函数曲线仅触及横轴和不穿越的零点无法确定函数曲线仅触及横轴和不穿越的零点 第二个输入量第二个输入量x0表示零点的初始猜测表示零点的初始猜测.他可以是标量或二元向他可以是标量或二元向量量.当当x0取标量时取标量时,该指令将在它两侧寻找一个与之最靠近的零该指令将在它两侧寻找一个与之最靠近的零点点;当取二元向量当取二元向量a,b时时,该指令将在区间该指令将在区间a,b内寻找一个零内寻找一
38、个零点点.输入宗量输入宗量options是优化迭代采用的参数是优化迭代采用的参数.P1,P2是向函数是向函数fun传传递的参数递的参数.z是所求零点的自变量值是所求零点的自变量值.f_z是函数值是函数值,exitflag表明程序终表明程序终止计算的条件止计算的条件.若若exitflag0,表明找到零点后退出表明找到零点后退出39例例4.4.1-1 求函数求函数f(x)=2x-10 x的零点的零点(1)绘制函数图形绘制函数图形x=-10:0.1:10;y=2.x-10*x;hold on;plot(x,y,r,LineWidth,5);plot(x,zeros(size(x),k);%便于查看零
39、点便于查看零点xlabel(x),ylabel(y);title(曲线方程式为：曲线方程式为：y=2x-10*x);(2)获取初步近似解获取初步近似解zoom on;tt,yy=ginput(2);zoom off;40func=2x-10*x;x01,err1,exitflag1=fzero(func,tt(1),);x02,err2,exitflag2=fzero(func,tt(2),);(3)利用利用fzero指令求精确解指令求精确解str=方程的根为方程的根为:x1=,num2str(x01),和和 x2=,num2str(x02);disp(str);414.4.2 求函数极值点求
40、函数极值点许多科学研究和工程计算问题都可以归结为一个极值问题许多科学研究和工程计算问题都可以归结为一个极值问题,如能量最小，时如能量最小，时间最短间最短,最佳拟合最佳拟合,最短路径等最短路径等.又因为又因为f(x)的极小值问题等价于的极小值问题等价于-f(x)的极大的极大值问题值问题,所以所以matlab只有处理极小值的指令只有处理极小值的指令确切地说确切地说,这里讨论的只是这里讨论的只是“局域极值局域极值”问题问题.“全域最小全域最小”问问题要复杂得多题要复杂得多.至今没有一个系统性的方法可求解一般的至今没有一个系统性的方法可求解一般的“全域全域最小最小”问题问题.对于一元二元函数对于一元二
41、元函数,作图观察对确定全域最小值有作图观察对确定全域最小值有很好的应用价值很好的应用价值.但更多元的函数但更多元的函数,就很难使用作图法就很难使用作图法.matlab求函数极值的三条指令求函数极值的三条指令x,fval,exitflag,output=fminbnd(fun,x1,x2,opt,p1,p2,)-求一元函数在区间求一元函数在区间(x1,x2)中极小值的指令的最完整格式中极小值的指令的最完整格式x,fval,exitflag,output=fminsearch(fun,x0,opt,p1,p2,)-单纯形法求多元函数极值点指令的最完整格式单纯形法求多元函数极值点指令的最完整格式x,
42、fval,exitflag,output,grad,hessian=fminunc(fun,x0,opt,p1,p2,)-拟牛顿法求多元函数极值点指令的最完整格式拟牛顿法求多元函数极值点指令的最完整格式42例例4.4.2-1 求例求例4.4.1-1 中的函数的最小值中的函数的最小值(1)函数采用字符串表达式函数采用字符串表达式func=2x-10*x;x1=0,x2=8;%通过观察通过观察,发现最小值在此区间发现最小值在此区间Xmin,Ymin=fminbnd(func,x1,x2)Xmin=Ymin=(2)函数采用函数采用M文件函数的句柄文件函数的句柄%M文件内容，保存为文件内容，保存为fu
43、nction y=f(x)y=2x-10*x;%M文件内容文件内容x1=0;x2=8;Hf=func;Xmin,Ymin=fminbnd(Hf,x1,x2)434.4.3 数值积分数值积分数值积分有闭型和开型算法之分两者的区别在于：是否需要计数值积分有闭型和开型算法之分两者的区别在于：是否需要计算积分区间端点处的函数值算积分区间端点处的函数值matlab中仅提供闭型数值积分指令中仅提供闭型数值积分指令:q=quad(fun,a,b,tol,trace,p1,p2,)-采用递推自采用递推自 适应适应Simpson法计算积分法计算积分q=qudal(fun,a,b,tol,trace,p1,p2,
44、)-采用递推自采用递推自适应适应Lobatlo法求数值积分法求数值积分SS=dblquad(fun,inmin,inmax,outmin,outmax,tol,method)-二重积分指令二重积分指令fun-被积函数被积函数,a,b-积分上下限积分上下限tol-绝对误差控制量，缺省时积分的绝对精度为绝对误差控制量，缺省时积分的绝对精度为1x10-6trace-取非时，将随积分的进程逐点画出被积函数取非时，将随积分的进程逐点画出被积函数44例求例求f(x)=3x2在区间在区间0,2的积分的积分y=3*x.2;a=0;b=2;q1=quad(y,a,b)q2=quadl(y,a,b)q1=q2=4
45、54.5 解常微分方程解常微分方程matlab为解常微分方程初值问题提供了一组配套齐全为解常微分方程初值问题提供了一组配套齐全,结构严结构严整的指令整的指令,包括包括:ode45,ode23,ode113,ode23t,ode15s,ode23s,ode23tb.在此只介绍最常用的在此只介绍最常用的ode45的基本使用方法的基本使用方法.ode45使用方法使用方法:x,y=ode45(file_name,xspan,y0)file_name=file_namex-微分区间的点系列微分区间的点系列y-原函数在微分区间点系列上的函数值原函数在微分区间点系列上的函数值file_name-档名档名,保
46、存的保存的M文件名文件名xspan-微分区间微分区间y0-初始值初始值46例例4.5-1 解微分方程：解微分方程：y=3x2,y(0)=1,微分区间微分区间0,5%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)474.6 数字滤波数字滤波y=filter(B,A,x,dim)-运用数字滤波对输入信号运用数字滤波对输入信号x滤波滤波滤波器的作用滤波器的作用:通过对采样数据信号进行数学运算处理来达到频域滤波的目的通过对采样
47、数据信号进行数学运算处理来达到频域滤波的目的.例例 设计一个低通滤波器设计一个低通滤波器,从受噪声干扰的多频率混合信号从受噪声干扰的多频率混合信号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;%Nyquest频率频率B,A=butter(10,30/wn);%截止频率为的截止频率为的10阶阶ButterWorth低通滤波器低通滤波器y=filter(B,A,x);%滤波处理滤波处理plot(t,x,b-,t,y,r.,MarkerSize,10);4849Ac=1,-0.5,-0.5;0,0.87,-0.87;Ap=cos(pi/4),sin(pi/4);-sin(pi/4),cos(pi/4);t=0:0.1:10;IA=220*cos(2*pi*60*t);IB=220*cos(2*pi*60*t-2*pi/3);IC=220*cos(2*pi*60*t+2*pi/3);I=IA;IB;IC;i=0.816*Ap*(Ac*I);P122-850