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

（三）Matlab的高级数值计算,关系运算逻辑运算多项式计算数值积分与微分数据插值曲线拟合方程组求解傅立叶分析,（三）Matlab的高级数值计算关系运算,matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。f(x)=anxn+an-1xn-1+a0 可用行向量 p=an an-1 a1+a0表示(1)poly 产生特征多项式系数向量特征多项式一定是n+1维的,3.多项式运算,matlab语言把多项式表达成一个行向量，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-72x-27的matlab描述方法，可用：p1=poly2str(p,x)函数文件，显示数学多项式的形式p1=x3-6 x2-72 x-27,例:a=1 2 3;4 5 6;7 8 0;,利用roots求多项式的根 r=roots(p)r=12.1229-5.7345-0.3884,当然可用poly令其返回多项式形式p2=poly(r)p2=1.00-6.00-72.00-27.00,matlab规定多项式系数向量用行向量表示，一组根用列向量表示。,(2)、多项式求根,利用roots求多项式的根当然可用poly令其返回多项式形式,求多项式零点:roots(p);,例1.求方程 x3 4x+5=0 的解P=1 0-4 5;roots(P)ans=-2.4567 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,求多项式零点:roots(p);例1.求方程 x3,(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-12089-43779,(3)Polyval 计算系数为p的多项式在标量或向量x处,(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,(4).conv多项式乘运算例:a(x)=x2+2x+3;,(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.00,conv(a,d)%因余数为零，可通过a,d的乘积返回 原多项式 ans=4 13 28 27 18,(5).deconv多项式除运算a=1 2 3;d,(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,(6).多项式微分matlab提供了polyder函数多项式,k=conv(p,q)k,r=deconv(p,q),k=polyder(p)k=polyder(p,q)k,d=polyder(p,q),y=polyval(p,x),x=roots(p),多项式运算小结,多项式运算中，使用的是多项式 系数向量，不涉及符号计算！,k=conv(p,q)k=polyd,4.数值微积分,4.1 数值积分 数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、欧拉法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等 基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,4.数值微积分4.1 数值积分,被积函数由一个表格定义 在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=10之间所围面积，并讨论步长和积分方法对精度的影响,建模：将x的被积区间分为n段，各段长度为xi（i=1，2，n+1）,欧拉法数值积分为 梯形法为,例：用数值积分法求 在x=0,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;%用梯形法求积分 stairs(x1,y1),plot(x1,y1)%画出欧拉法及梯形法的积分区域 dx,s,q,pause(1),hold off%给出不同步长下的结果 end,MATLAB程序,运行结果：步长 欧拉法解 梯形法解 2 910 810 1 865 815 0.5 841.25 816.25 0.1 821.65 816.65,解析法得到的精确解为2450/3=816.6667,右图显示曲线的积分面积，在曲线斜率为负的情况下，欧拉法解偏大，梯形法偏小，精确解位于两者之间,步长相同时，梯形法的精度比欧拉法要高,运行结果：解析法得到的精确解为2450/3=816.6667,变步长辛普生法 基于变步长辛普生法，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 function y=ex5f(x)y=-x.*x+115;应用quad函数 s=quad（ex5f,0,10）运行结果：s=816.6667,变步长辛普生法上例中，建立函数文件ex5f.m,牛顿柯特斯法 基于牛顿柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：I,n=quad8(fname,a,b,tol,trace)其中参数的含义和quad函数相似，只是tol的缺省值取10-6。该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。,牛顿柯特斯法,例：分别用quad函数和quad8函数求定积分的 近似值，并在相同的积分精度下，比较函数 的调用次数。调用函数quad求定积分：format long;fx=inline(exp(-x);%定义函数I,n=quad(fx,1,2.5,1e-10)结果：I=0.28579444254766n=65,例：分别用quad函数和quad8函数求定积分的,调用函数quad8求定积分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)运行结果I=0.28579444254754n=33,调用函数quad8求定积分：,二重定积分的数值求解 使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。,二重定积分的数值求解,例 计算二重定积分(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.57449318974494,例 计算二重定积分,4.2 数值微分 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，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，按行计算差分。,4.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 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 0,例 生成以向量V=1,2,3,4,5,6为基础的范得蒙矩,5.数据插值,5.数据插值,插值方法,一维插值的定义已知n个节点，求任意点处的函数值。分段线性插值 linear多项式插值 cubic样条插值 spline最近邻插值方法 nearesty=interp1(x0,y0,x,method),插值方法一维插值的定义已知n个节点，求任意点处的函数值。,举例 多项式插值原理,已知 f(x)在点 xi 上的函数值 yi=f(xi),(i=0,1,2,n),求多项式 P(x)=a0+a1x+anxn满足:P(xk)=yk(k=0,1,n),举例 多项式插值原理已知 f(x)在点 xi 上的函数值,一维插值:yi=interp1(x,y,xi,method),例：x=0:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);%默认线性插值plot(x,y,o,xi,yi),一维插值:yi=interp1(x,y,xi,第3讲Matlab的高级计算310课件,%已知数据t=1900:10:1990;p=75.995 91.972 105.711 123.203 131.669.150.697 179.323 203.212 226.505 249.633;%使用不同的方法进行插值运算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 Spline);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);,对某城市人口数量统计统计插值,%已知数据对某城市人口数量统计统计插值,第3讲Matlab的高级计算310课件,二维插值高维插值的一种，主要应用于图像处理和数据的可视化，对两个变量的函数z=f(x,y)进行插值,zi=interp2(x,y,z,xi,yi,method)method 用法类似一维插值，只是参数linear表示的是双线性插值方法,二维插值zi=interp2(x,y,z,xi,yi,例：以二元函数 为例，首先获取基础数据，而后使用二维插值得到更细致的数据，完成绘图,%生成基础测量数据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,%进行二维插值运算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向量数据变为矩阵数据，生成网格点,%生成基础测量数据%进行二维插值运算,其他插值方法：除了MATLAB内置的插值函数以外，还可以根据工程实践的需要，编写适用于具体领域的插值函数文件来加以调用，如：牛顿插值方法 拉格朗日（Lagrange）插值 Chebyshev多项式插值方法,其他插值方法：,6.多项式拟和,在科学与工程领域，曲线拟合的主要功能是寻求平滑的曲线来最好的表现带有噪声的测量数据，从测量数据中找到变量函数的变化趋势，最后得到曲线拟合的函数表达式,6.多项式拟和在科学与工程领域，曲线拟合的主要功能是寻求平,已知数据表,x x1 x2 xmf(x)y1 y2 ym,求拟合函数:(x)=a0+a1x+anxn使得,离散数据的多项式拟合,达到最小,已知数据表 x x1 x2,设 i=1,2,N表示按拟合直线 求得的近似值，一般地说，它不同于实测值，两者之差称 残差。显然，残差的大小是衡量拟合好坏的重要标志，具体地说，可以采用下列三种准则：,使残差的最大绝对值为最小：使残差的绝对值之和最小：使残差的平方和为最小：,直线拟合,设 i=1,2,对于给顶的数据点,直线拟合问题可用数学语言描述如下：,求作一次式,使总误差,最小。,有时候所给出数据点用直线拟合不合适，这时可考虑用多项式拟合，用数学语言描述如下：对于给定的一组数据，寻求作次多项式（）使总误差 为最小。,曲线拟合,有时候所给出数据点用直线拟合不合适，这时可考虑用多项式拟合，,多项式拟合命令,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处的值,多项式拟合命令P=polyfit(X,Y,n)求出(最小,例：如何预报人口的增长 人口的增长是当前世界上引起普遍关注的问题，并且我们会发现在不同的刊物预报同一时间的人口数字不相同，这显然是由于用了不同的人口模型计算的结果。,例如：1949年1994年我国人口数据资料如下：年 份 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=0.146,介绍两个简单模型模,模型二：指数增长模型,（用简单的线性最小二乘法）,用Matlab软件计算得：a=2.33,b=0.0179即:,模型二：指数增长模型（用简单的线性最小二乘,程序如下：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=exp(b(2)*exp(b(1)*x1);plot(x,y,*)hold onplot(x1,y1,-r)hold on plot(x1,y2,-k)legend(原曲线,模型一曲线,模型二曲线),程序如下：x=1949 1954 1959 1964 1,程序执行后得到下面图形,程序执行后得到下面图形,结论的比较如下表：,结论的比较如下表：年 份 xi194919541959196,结果分析：,模型 I 2019年13.43亿，2019年14.16亿,误差Q1=0.2915 0.7437=Q2.线性模型更适合中 国人口的增长。2.预报：2019年12.55亿，13.43亿3.人口白皮书：2019年13.3亿，2019年14亿,模型II 2019年14.94亿，2019年16.33亿,结果分析：模型 I 2019年13.43,曲线拟合图形界面,例：在命令窗口输入拟合的基础数据x=1 3 5 7 9 2 4 6 8 10;y=0.0881 0.929 2.4932 4.9292 7.9605.0.9536 2.4836 3.4173 6.3903 10.2443;plot(x,y,ro)grid on在图形窗口tools中选择Basic fitting,曲线拟合图形界面例：,7.方程组求解,7.方程组求解,对于方程ax+b，a 为anm矩阵，有三种情况：当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程 matlab定义的除运算可以很方便地解上述三种方程,通过除运算左除和右除来解方程组,对于方程ax+b，a 为anm矩阵，有三种情通过除运算左除,恰定方程组的解,方程ax+b(a为非奇异)x=a-1 b 矩阵逆两种解:x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程,恰定方程组的解方程ax+b(a为非奇异),方程ax=ba=1 2;2 3;b=8;13;x=inv(a)*b x=ab x=x=2.00 2.00 3.00 3.00,=,a x=b,例:x1+2x2=8 2x1+3x2=13,方程ax=b,线性方程组求解,线性方程组求解,linsolve(A,b)：解线性方程组,例：解方程组,A=1 2 1;1 0 1;1 3 0;b=2;3;8;x=linsolve(A,b),b是列向量！,线性方程组求解 线性方程组求解linsolve(A,b)：解,非线性方程的根,Matlab 非线性方程的数值求解,fzero(f,x0)：求方程 f=0 在 x0 附近的根。,方程可能有多个根，但 fzero 只给出距离 x0 最近的一个,fzero 先找出一个包含 x0 的区间，使得 f 在这个区间两个端点上的函数值异号，然后再在这个区间内寻找方程 f=0 的根；如果找不到这样的区间，则返回 NaN。,x0 是一个标量，不能缺省,非线性方程的根 Matlab 非线性方程的数值求解fzero,非线性方程的根,fzero 的另外一种调用方式,fzero(f,a,b),方程在 a,b 内可能有多个根，但 fzero 只给出一个,求方程 f=0 在 a,b 区间内的根。,参数 f 可通过以下方式给出：,fzero(x3-3*x+1,2);f=inline(x3-3*x+1);fzero(f,2)fzero(x)x3-3*x+1,2);,f 不是方程！也不能使用符号表达式！,非线性方程的根 fzero 的另外一种调用方式fzero(f,例：,fzero(sin(x),10)fzero(sin,10),fzero(x3-3*x+1,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 求零点时可以先通过作图确定零点的大致范围,例：fzero(sin(x),10)fzer,8.傅立叶分析,傅立叶变换是分析周期或非周期信号的频率特性的数学工具，在信号处理领域有着广泛应用,8.傅立叶分析傅立叶变换是分析周期或非周期信号的频率特性的数,例：分析周期信号 其中w1=50HZ,w2=120HZ,t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(size(t);plot(1000*t(1:50),y(1:50),例：分析周期信号,%进行傅立叶变换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),N=find(Pyy20);%找出频率较弱的信号,