第3讲Matlab的高级计算310课件.ppt

1、（三）Matlab的高级数值计算 关系运算关系运算 逻辑运算逻辑运算 多项式计算多项式计算 数值积分与微分数值积分与微分 数据插值数据插值 曲线拟合曲线拟合 方程组求解方程组求解 傅立叶分析傅立叶分析matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。f(x)=anxn+an-1xn-1+a0 可用行向量 p=an an-1 a1+a0表示(1)poly 产生特征多项式系数向量特征多项式一定是n+1维的3.多项式运算 例:a=1 2 3;4 5 6;7 8 0;p=poly(a)p=1.00 -6.00 -72.00 -27.00这是多项式p(x)=x3-6x2-72

2、x-27的matlab描述方法，可用：p1=poly2str(p,x)函数文件，显示数学多项式的形式p1=x3-6 x2-72 x-27利用roots求多项式的根 r=roots(p)r=12.1229 -5.7345 -0.3884当然可用当然可用poly令其返回多项式形式令其返回多项式形式p2=poly(r)p2=1.00 -6.00 -72.00 -27.00vmatlab规定多项式系数向量用行向量表示，一组根用列向量表示。(2)、多项式求根 求多项式零点求多项式零点:roots(p);例例1.求方程求方程 x3 4x+5=0 的解的解P=1 0-4 5;roots(P)ans=-2.4

3、567 1.2283+0.7256i 1.2283-0.7256i例例2 求求 x3 8x2+6x 30=0的解的解P=1-8 6-30;r=roots(P)r=7.7260 0.1370+1.9658i 0.1370-1.9658i(3)Polyval 计算系数为p的多项式在标量或向量x处的值 X=pascal(4)X=1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 p=poly(X)p=1 -29 72 -29 1polyval(p,X)ans=16 16 16 16 16 15 -140 -563 16 -140 -2549 -12089 16 -563 -1208

4、9 -43779(4).conv多项式乘运算例:a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=1 2 3;b=4 5 6;c=conv(a,b)c=4.00 13.00 28.00 27.00 18.00p=poly2str(c,x)p=4 x4+13 x3+28 x2+27 x+18(5).deconv多项式除运算a=1 2 3;c=4.00 13.00 28.00 27.00 18.00d=deconv(c,a)d=4.00 5.00 6.00d,r=deconv(c,a)余数余数c除以除以a后的商后的商conv(a,d)%因余数为零，

5、可通过a,d的乘积返回 原多项式 ans=4 13 28 27 18(6).多项式微分matlab提供了polyder函数多项式的微分。命令格式：polyder(p):求p的微分polyder(a,b):求多项式a,b乘积的微分p,q=polyder(a,b):求多项式a,b商的微分例：a=1 2 3 4 5;poly2str(a,x)ans=x4+2 x3+3 x2+4 x+5b=polyder(a)b=4 6 6 4poly2str(b,x)ans=4 x3+6 x2+6 x+4 k=conv(p,q)k,r=deconv(p,q)k=polyder(p)k=polyder(p,q)k,d

6、=polyder(p,q)y=polyval(p,x)x=roots(p)多项式运算小结多项式运算中，多项式运算中，使用的是多项式使用的是多项式 系数向量系数向量，不涉及符号计算！不涉及符号计算！4.数值微积分数值微积分4.1 数值积分数值积分 数值积分基本原理数值积分基本原理 求解定积分的数值方法多种多样，如简单的求解定积分的数值方法多种多样，如简单的梯梯形法形法、欧拉法、辛普生欧拉法、辛普生(Simpson)法法、牛顿牛顿柯特斯柯特斯(Newton-Cotes)法等法等 基本思想都是将整个积分区间基本思想都是将整个积分区间a,b分成分成n个子个子区间区间xi,xi+1，i=1,2,n，其中

7、，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。这样求定积分问题就分解为求和问题。被积函数由一个表格定义被积函数由一个表格定义 在在MATLAB中，对由表格形式定义的函数关系中，对由表格形式定义的函数关系的求定积分问题用的求定积分问题用trapz(X,Y)函数。其中向量函数。其中向量X,Y定义函数关系定义函数关系Y=f(X)。例例 用用trapz函数计算定积分。函数计算定积分。命令如下：命令如下：X=1:0.01:2.5;Y=exp(-X);%生成函数关系数据向量生成函数关系数据向量trapz(X,Y)ans=0.28579682416393 例：用数值积分法求 在x=0到x=

8、10之间所围面积，并讨论步长和积分方法对精度的影响 115xy2建模：将x的被积区间分为n段，各段长度为xi（i=1，2，n+1）,欧拉法数值积分为 梯形法为 n1iix)i(ysin1ix2)1i(y)i(y(q MATLAB程序 for dx=2,1,0.5,0.1%设不同步长 x=0:.1:10;y=-x.*x+115;plot(x,y,g),hold%画出被积曲线 x1=0:dx:10;%根据不同步长取样点 y1=-x1.*x1+115;%求取样点上的y1 n=length(x1);s=sum(y1(1:n-1)*dx;%用欧拉法求积分 q=trapz(y1)*dx;%用梯形法求积分

9、stairs(x1,y1),plot(x1,y1)%画出欧拉法及梯形法的积分区域 dx,s,q,pause(1),hold off%给出不同步长下的结果 end 运行结果：步长 欧拉法解 梯形法解 2 910 810 1 865 815 0.5 841.25 816.25 0.1 821.65 816.65解析法得到的精确解为2450/3=816.6667012345678910020406080100120右图显示曲线的积分面积，在曲线斜率为负的情况下，欧拉法解偏大，梯形法偏小，精确解位于两者之间 步长相同时，梯形法的精度比欧拉法要高变步长辛普生法变步长辛普生法 基于变步长辛普生法，基于变步

10、长辛普生法，MATLAB给出了给出了quad函数来求定积分。该函数来求定积分。该函数的调用格式为：函数的调用格式为：I,n=quad(fname,a,b,tol,trace)其中其中fname是被积函数名。是被积函数名。a和和b分别是定积分的下限和上限。分别是定积分的下限和上限。tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分控制是否展现积分过程，若取非过程，若取非0则展现积分过程，取则展现积分过程，取0则不展现，缺省时取则不展现，缺省时取trace=0。需要事先建立函数文件需要事先建立函数文件fname上例中，建立函数文件ex5f.m

11、function y=ex5f(x)y=-x.*x+115;应用quad函数 s=quad（ex5f,0,10）运行结果：s=816.6667 牛顿柯特斯法牛顿柯特斯法 基于牛顿柯特斯法，基于牛顿柯特斯法，MATLAB给出了给出了quad8函数来求定积分。该函数的调用格函数来求定积分。该函数的调用格式为：式为：I,n=quad8(fname,a,b,tol,trace)其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。该函数可以更精确地该函数可以更精确地求出定积分的值，且一般情况下函数调求出定积分的值，且一般情况下函数调用的步数明显小于用的

12、步数明显小于quad函数，从而保证函数，从而保证能以更高的效率求出所需的定积分值。能以更高的效率求出所需的定积分值。例例：分别用分别用quad函数和函数和quad8函数求定积分的函数求定积分的 近似值，并在相同的积分精度下，比较函数近似值，并在相同的积分精度下，比较函数 的调用次数。的调用次数。调用函数调用函数quad求定积分：求定积分：format long;fx=inline(exp(-x);%定义函数定义函数I,n=quad(fx,1,2.5,1e-10)结果：结果：I=0.28579444254766n=65调用函数调用函数quad8求定积分：求定积分：format long;fx=i

13、nline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)运行结果运行结果I=0.28579444254754n=33二重定积分的数值求解二重定积分的数值求解 使用使用MATLAB提供的提供的dblquad函数就可函数就可以直接求出上述二重定积分的数值解。以直接求出上述二重定积分的数值解。该该函数的调用格式为：函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求该函数求f(x,y)在在a,bc,d区域上的二重定积区域上的二重定积分。分。参数参数tol，trace的用法与函数的用法与函数quad完全相同。完全相同。例例 计算二重定积分计算二重

14、定积分(1)建立一个函数文件建立一个函数文件fxy.m：function f=fxy(x,y)f=exp(-x.2/2).*sin(x.2+y);(2)调用调用dblquad函数求解。函数求解。I=dblquad(fxy,-2,2,-1,1)I=1.57449318974494112222x)yx(sine24.2 数值微分数值微分 在在MATLAB中，没有直接提供求数值导中，没有直接提供求数值导数的函数，只有计算向前差分的函数数的函数，只有计算向前差分的函数diff，其调用格式为：其调用格式为：DX=diff(X)：计算向量：计算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i

15、)，i=1,2,n-1。DX=diff(X,n)：计算：计算X的的n阶向前差分。例阶向前差分。例如，如，diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)：计算矩阵：计算矩阵A的的n阶差分，阶差分，dim=1时时(缺省状态缺省状态)，按列计算差分；，按列计算差分；dim=2，按行计算差分。，按行计算差分。例例 生成以向量生成以向量V=1,2,3,4,5,6为基础的范得蒙矩为基础的范得蒙矩阵，按列进行差分运算。阵，按列进行差分运算。命令如下：命令如下：V=vander(1:6)DV=diff(V)%计算计算V的一阶差分的一阶差分V=1 1 1 1 1 1 32 16

16、8 4 2 1 243 81 27 9 3 1 1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1DV=31 15 7 3 1 0 211 65 19 5 1 0 781 175 37 7 1 0 2101 369 61 9 1 0 4651 671 91 11 1 05.数据插值插值方法一维插值的定义已知n个节点，求任意点处的函数值。分段线性插值 linear多项式插值 cubic样条插值 spline最近邻插值方法 nearesty=interp1(x0,y0,x,method)举例举例 多项式插值原理多项式插值原理已知已

17、知 f(x)在点在点 xi 上的函数值上的函数值 yi=f(xi),(i=0,1,2,n)求多项式求多项式 P(x)=a0+a1x+anxn满足满足:P(xk)=yk (k=0,1,n)P10(t)f(t)f(x)一维插值一维插值:yi=interp1(x,y,xi,method)methodnearest 最近邻点插值最近邻点插值 linear 线性插值线性插值 spline 样条插值样条插值 cubic 三次多项式三次多项式 插值插值例：例：x=0:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);%默认线性插值默认线性插值plot(x,y,o,xi,y

18、i)012345678910-1-0.8-0.6-0.4-0.200.20.40.60.81linear012345678910-1-0.8-0.6-0.4-0.200.20.40.60.81cubic012345678910-1-0.8-0.6-0.4-0.200.20.40.60.81nearest012345678910-1-0.8-0.6-0.4-0.200.20.40.60.81spline%已知数据t=1900:10:1990;p=75.995 91.972 105.711 123.203 131.669.150.697 179.323 203.212 226.505 249.63

19、3;%使用不同的方法进行插值运算x=1900:0.01:1990;yi_linear=interp1(t,p,x);yi_spline=interp1(t,p,x,spline);yi_cubic=interp1(t,p,x,cubic);yi_nearest=interp1(t,p,x,nearest);%绘制对比图形subplot(2,1,1);plot(t,p,ko);hold on;plot(x,yi_linear,g,LineWidth,1.5);hold onplot(x,yi_spline,y,LineWidth,1.5);grid on;title(Linear Vs Spli

20、ne);subplot(2,1,2);plot(t,p,ko);hold on;plot(x,yi_cubic,m,LineWidth,1.5);hold onplot(x,yi_nearest,k,LineWidth,1);grid on;title(Cubic Vs nearest);对某城市人口数量统计统计插值190019101920193019401950196019701980199050100150200250Linear Vs Spline190019101920193019401950196019701980199050100150200250Cubic Vs nearest二

21、维插值二维插值高维插值的一种，主要应用于图像处理和数据的可视化，对两个变量的函数z=f(x,y)进行插值zi=interp2(x,y,z,xi,yi,method)method 用法类似一维插值，只是参数用法类似一维插值，只是参数linear表示的是双线性插值方法表示的是双线性插值方法 例：以二元函数 为例，首先获取基础数据，而后使用二维插值得到更细致的数据，完成绘图2222yxyxsin)y,x(z%生成基础测量数据x=-3*pi:3*pi;y=x;X,Y=meshgrid(x,y);R=sqrt(X.2+Y.2)+eps;Z=sin(R)./R;%eps正的极小值，防止出 现0/0-10-

22、50510-10-50510-0.500.510.10.20.30.40.50.60.70.80.9%进行二维插值运算xi=linspace(-3*pi,3*pi,100);yi=linspace(-3*pi,3*pi,100);XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,cubic);%meshgrid向量数据变为矩阵数据，生成网格点-10-50510-10-50510-0.500.51-0.200.20.40.60.8 其他插值方法：除了MATLAB内置的插值函数以外，还可以根据工程实践的需要，编写适用于具体领域的插值函数文件来加以调用，如：牛

23、顿插值方法 拉格朗日（Lagrange）插值 Chebyshev多项式插值方法6.多项式拟和在科学与工程领域，曲线拟合的主要功能是寻求平滑的曲线来最好的表现带有噪声的测量数据，从测量数据中找到变量函数的变化趋势，最后得到曲线拟合的函数表达式已知数据表已知数据表 x x1 x2 xmf(x)y1 y2 ym求拟合函数求拟合函数:(x)=a0+a1x+anxn使得使得离散数据的多项式拟合离散数据的多项式拟合mkknknkyxaxaa1210)(达到最小达到最小设 i=1,2,N表示按拟合直线 求得的近似值，一般地说，它不同于实测值，两者之差称 残差。显然，残差的大小是衡量拟合好坏的重要标志，具体地

24、说，可以采用下列三种准则：使残差的最大绝对值为最小：使残差的绝对值之和最小：使残差的平方和为最小：maxminiie miniie 2miniieiiyabxyabxiiieyy直线拟合直线拟合 (xi,yi)1,2,3,iNyabx1()NiiiQyabx对于给顶的数据点直线拟合问题可用数学语言描述如下：,求作一次式,使总误差最小。有时候所给出数据点用直线拟合不合适，这时可考虑用多项式拟合，用数学语言描述如下：对于给定的一组数据 ，寻求作次多项式（）使总误差 为最小。(,)iix y1,2,3,iNmNmjjj oya x21()Nmjjiij oQya x曲线拟合曲线拟合多项式拟合命令多项

25、式拟合命令P=polyfit(X,Y,n)求出求出(最小二乘意义下最小二乘意义下)n次拟合多项式次拟合多项式P(x)=a0 xn+a1xn-1+an-1x+an计算结果为系数计算结果为系数:P=a0,a1,an-1,an 多项式求值命令多项式求值命令:y1=polyval(P,x)其中其中,P是是n次多项式的系数次多项式的系数,x是自变量的值是自变量的值,y1是多是多项式在项式在x处的值处的值例：如何预报人口的增长例：如何预报人口的增长 人口的增长是当前世界上引起普遍关注的问题，并且我们会发现在不同的刊物预报同一时间的人口数字不相同，这显然是由于用了不同的人口模型计算的结果。例如：1949年1

26、994年我国人口数据资料如下：年 份 xi 1949 1954 1959 1964 1969 1974 1979 19841989 1994人口数 yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 建模分析我国人口增长的规律,预报2019年我国人口数。介绍两个简单模型介绍两个简单模型模型一：假设：人口随时间线性地增加模型一：假设：人口随时间线性地增加拟合的精度:Q=ei 2=(yi-a b xi)2,误差平方和。模型：y=1.93+0.146 x参数估计 观测值的模型：yi=a+b xi+ei，i=1,n模型：y=a+b x可以算出：a=1.93，b=

27、0.146 模型二：模型二：指数增长模型指数增长模型（用简单的线性最小二乘法）（用简单的线性最小二乘法）用用Matlab软件计算得：软件计算得：a=2.33,b=0.0179即即:0.01792.33xyebxaey bxlnalny 程序如下：程序如下：x=1949 1954 1959 1964 1969 1974 1979 1984 1989 1994;y=5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8;a=polyfit(x,y,1);x1=1949:10:1994;y1=a(2)+a(1)*x1;b=polyfit(x,log(y),1);y2=e

28、xp(b(2)*exp(b(1)*x1);plot(x,y,*)hold onplot(x1,y1,-r)hold on plot(x1,y2,-k)legend(原曲线,模型一曲线,模型二曲线)程序执行后得到下面图形程序执行后得到下面图形结论的比较如下表：结论的比较如下表：年 份 xi1949195419591964196919741979198419891994人口数 yi5.466.778.19.19.810.311.311.8模型一值5.245.976.77.438.168.99.6210.3611.0911.82误 差0.160.030-0.43-0.060.20.18-0.060.

29、01-0.02模型模型二值二值5.555.556.066.066.626.627.237.237.97.98.648.649.449.4410.3110.3111.2611.2612.3112.31误差误差-0.15-0.060.08-0.230.20.460.36-0.01-0.13-0.51结果分析：结果分析：模型模型 I 2019年年13.43亿，亿，2019年年14.16亿亿1.误差误差Q1=0.2915 m时，此方程成为“超定”方程 当n A=1 2 1;1 0 1;1 3 0;b=2;3;8;x=linsolve(A,b)22338xyzxzxy b是列向量！是列向量！非线性方程的

30、根q Matlab 非线性方程的数值求解非线性方程的数值求解fzero(f,x0)：求方程求方程 f=0 在在 x0 附近的根。附近的根。l 方程可能有多个根，但方程可能有多个根，但 fzero 只给出距离只给出距离 x0 最近的一个最近的一个l fzero 先找出一个包含先找出一个包含 x0 的区间，使得的区间，使得 f 在这个区间在这个区间两个端点上的函数值异号，然后再在这个区间内寻找方程两个端点上的函数值异号，然后再在这个区间内寻找方程 f=0 的根；如果找不到这样的区间，则返回的根；如果找不到这样的区间，则返回 NaN。l x0 是一个标量，不能缺省是一个标量，不能缺省非线性方程的根q

31、 fzero 的另外一种调用方式的另外一种调用方式fzero(f,a,b)l 方程在方程在 a,b 内可能有多个根，但内可能有多个根，但 fzero 只给出一个只给出一个l 求方程求方程 f=0 在在 a,b 区间内区间内的根。的根。q 参数参数 f 可通过以下方式给出：可通过以下方式给出：l fzero(x3-3*x+1,2);l f=inline(x3-3*x+1);fzero(f,2)l fzero(x)x3-3*x+1,2);l f 不是方程！也不能使用符号表达式！不是方程！也不能使用符号表达式！例：fzero(sin(x),10)fzero(sin,10)fzero(x3-3*x+1

32、,1)fzero(x3-3*x+1,1,2)fzero(x3-3*x+1=0,1)X fzero(x3-3*x+1,-2,0)f=inline(x3-3*x+1);fzero(f,-2,0)用用 fzero 求零点时可以先通过作图确定零点的大致范围求零点时可以先通过作图确定零点的大致范围8.傅立叶分析 傅立叶变换是分析周期或非周期信号的频率特性的数学工具，在信号处理领域有着广泛应用05101520253035404550-5-4-3-2-1012345)tw2(sin)tw2(sinx21例：分析周期信号 其中w1=50HZ,w2=120HZt=0:0.001:0.6;x=sin(2*pi*5

33、0*t)+sin(2*pi*120*t);y=x+2*randn(size(t);plot(1000*t(1:50),y(1:50)05010015020025030035040045050001020304050607080%进行傅立叶变换Y=fft(y,512);%512点的fftPyy=Y.*conj(Y)/512;%conj复共轭函数f=1000*(0:256)/512;%绘图figure,plot(f,Pyy(1:257)N=find(Pyy20);%找出频率较弱的信号 Pyy(N)=0;%去除 figure,plot(f,Pyy(1:257)05010015020025030035040045050001020304050607080