第三讲 Matlab的高级计算3.10.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-72x-27的matlab描述方法， 可用： p1=p

2、oly2str(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.00,matlab规定多项式系数向量用行向量表示，一组根用列向量表示。,(2) 、多项式求根,求多项式零点: roots(p);,例1.求方程 x3 4x + 5 = 0 的解 P=1 0 -4 5; roots(P) ans = -2.4567 1.2283 + 0.72

3、56i 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 -

4、563-12089-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.00 p=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.00 d=deconv(c,a) d =4.00 5.00 6.00,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 + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4,k = conv(

6、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),多项式运算小结,多项式运算中， 使用的是多项式 系数向量， 不涉及符号计算！,4.数值微积分,4.1 数值积分 数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、欧拉法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等 基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问

7、题。,被积函数由一个表格定义 在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）, 欧拉法数值积分为 梯形法为,MATLAB程序 for dx=2,1,0.5,0.1 % 设不同

8、步长 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,运行结果： 步长 欧拉法解 梯形法解 2 910 810 1 865 815 0.5

9、 841.25 816.25 0.1 821.65 816.65,解析法得到的精确解为2450/3=816.6667,右图显示曲线的积分面积，在曲线斜率为负的情况下，欧拉法解偏大，梯形法偏小，精确解位于两者之间,步长相同时，梯形法的精度比欧拉法要高,变步长辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为： I,n=quad(fname,a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。 tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时

10、取trace=0。 需要事先建立函数文件fname,上例中，建立函数文件ex5f.m 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。该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。,例： 分别

11、用quad函数和quad8函数求定积分的 近似值，并在相同的积分精度下，比较函数 的调用次数。 调用函数quad求定积分： format long; fx=inline(exp(-x); %定义函数 I,n=quad(fx,1,2.5,1e-10) 结果： I = 0.28579444254766 n = 65,调用函数quad8求定积分： format long; fx=inline(exp(-x); I,n=quad8(fx,1,2.5,1e-10) 运行结果 I = 0.28579444254754 n = 33,二重定积分的数值求解 使用MATLAB提供的dblquad函数就可以直接求

12、出上述二重定积分的数值解。 该函数的调用格式为： 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

13、， 其调用格式为： 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，按行计算差分。,例 生成以向量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

14、1024 256 64 16 4 1 3125 625 125 25 5 1 7776 1296 216 36 6 1 DV = 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,5.数据插值,插值方法,一维插值的定义已知n个节点，求任意点处的函数值。 分段线性插值 linear 多项式插值 cubic 样条插值 spline 最近邻插值方法 nearest y=interp1(x0,y0,x,method),举例 多项式插值原理,已知 f(x) 在点 xi 上的函数值 yi=

15、f(xi), (i=0,1,2,n),求多项式 P(x)=a0 + a1x + anxn 满足: P(xk)= yk (k = 0,1,n),一维插值: 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),%已知数据 t = 1900:10:1990; p = 75.995 91.972 105.711 123.203 131.669. 150.697 179.323 203.212 226.505 249.6

16、33; %使用不同的方法进行插值运算 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 on plot(x,yi_spline,y,LineWidth,1.5); grid on; titl

17、e(Linear Vs Spline); subplot(2,1,2); plot(t,p,ko);hold on; plot(x,yi_cubic,m,LineWidth,1.5);hold on plot(x,yi_nearest,k,LineWidth,1); grid on; title(Cubic Vs nearest);,对某城市人口数量统计统计插值,二维插值 高维插值的一种，主要应用于图像处理和数据的可视化，对两个变量的函数z=f(x,y)进行插值,zi=interp2(x, y, z, xi, yi, method) method 用法类似一维插值，只是参数linear表示的是

18、双线性插值方法,例：以二元函数 为例，首先获取基础数据，而后使用二维插值得到更细致的数据，完成绘图,%生成基础测量数据 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向量数据变为矩阵数据，生成网格点,

19、其他插值方法： 除了MATLAB内置的插值函数以外，还可以根据工程实践的需要，编写适用于具体领域的插值函数文件来加以调用， 如： 牛顿插值方法 拉格朗日（Lagrange）插值 Chebyshev多项式插值方法,6. 多项式拟和,在科学与工程领域，曲线拟合的主要功能是寻求平滑的曲线来最好的表现带有噪声的测量数据，从测量数据中找到变量函数的变化趋势，最后得到曲线拟合的函数表达式,已知数据表,求拟合函数: (x) = a0 + a1x + + anxn 使得,离散数据的多项式拟合,达到最小,设 i=1,2,N表示按拟合直线 求得的近似值，一般地说，它不同于实测值，两者之差称 残差。显然，残差的大小

20、是衡量拟合好坏的重要标志，具体地说，可以采用下列三种准则：,使残差的最大绝对值为最小： 使残差的绝对值之和最小： 使残差的平方和为最小：,直线拟合,对于给顶的数据点,直线拟合问题可用数学语言描述如下：,求作一次式,使总误差,最小。,有时候所给出数据点用直线拟合不合适，这时可考虑用多项式拟合，用数学语言描述如下： 对于给定的一组数据 ， ，寻求作次多项式（ ）使总误差 为最小。,曲线拟合,多项式拟合命令,P=polyfit(X,Y,n),求出(最小二乘意义下)n次拟合多项式 P(x)=a0 xn + a1xn-1 + + an-1x + an 计算结果为系数: P= a0, a1, , an-1

21、, an ,多项式求值命令: y1=polyval(P, x) 其中,P是n次多项式的系数,x是自变量的值,y1是多项式在x处的值,例：如何预报人口的增长 人口的增长是当前世界上引起普遍关注的问题，并且我们会发现在不同的刊物预报同一时间的人口数字不相同，这显然是由于用了不同的人口模型计算的结果。,例如：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 建模分析我国人口增长的规律,预报1999

22、年我国人口数。,介绍两个简单模型,模型一：假设：人口随时间线性地增加,拟合的精度: 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

23、.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(原曲线,模型一曲线,模型二曲线),程序执行后得到下面图形,结论的比较如下表：,结果分析：,模型 I 2005年13.43亿， 2010年14.16亿,误差Q1 = 0.2915 0.7437 = Q2.

24、 线性模型更适合中 国人口的增长。 2. 预报：1999年12.55亿，13.43亿 3. 人口白皮书： 2005年13.3亿， 2010年14亿,模型II 2005年14.94亿， 2010年16.33亿,曲线拟合图形界面,例： 在命令窗口输入拟合的基础数据 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. 方程组求解,对于方程ax+b，a 为an

25、m矩阵，有三种情 况： 当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程 matlab定义的除运算可以很方便地解上 述三种方程,通过除运算左除和右除来解方程组,恰定方程组的解,方程ax+b(a为非奇异) x=a-1 b 矩阵逆 两种解: x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程,方程ax=b a=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,线性方程组求解,线性方程组求解,linsolve(A,b)：解线性方程组,例：解方程组, A=1 2 1; 1 0 1; 1 3 0; b=2;3;8; x=linsolve(A,b),b是列向量！,非线性方程的根,Matlab 非线性方程的数值求解,fzero(f,x0)：求方程 f=0 在 x0 附近的根。,方程可能有多个根，