《Matlab的高级计算课件》由会员分享，可在线阅读，更多相关《Matlab的高级计算课件（59页珍藏版）》请在金锄头文库上搜索。
1、（三）Matlab的高级数值计算关系运算关系运算逻辑运算逻辑运算多项式计算多项式计算数值积分与微分数值积分与微分数据插值数据插值曲线拟合曲线拟合方程组求解方程组求解傅立叶分析傅立叶分析Matlab的高级计算matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。f(x)=anxn+an-1xn-1+a0可用行向量p=anan-1a1+a0表示(1)poly产生特征多项式系数向量特征多项式一定是n+1维的3.多项式运算Matlab的高级计算例:a=123;456;780;p=poly(a)p=1.00-6.00-72.00-27.00这是多项式p(x)=x3-6x2-72x
2、-27的matlab描述方法，可用：p1=poly2str(p,x)函数文件，显示数学多项式的形式p1=x3-6x2-72x-27Matlab的高级计算利用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) 、多项式求根 Matlab的高级计算求多项式零点求多项式零点: roots(p);例例1.求方程求方程 x3 4x + 5 = 0 的解的解P=
3、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.9658iMatlab的高级计算(3)Polyval计算系数为p的多项式在标量或向量x处的值X=pascal(4)X=1111123413610141020p=poly(X)p=1-2972-291polyval(p,X)ans=161616161615-140-56316-140-2549-1
4、208916-563-12089-43779Matlab的高级计算(4).conv多项式乘运算例:a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=123;b=456;c=conv(a,b)c=4.0013.0028.0027.0018.00p=poly2str(c,x)p=4x4+13x3+28x2+27x+18Matlab的高级计算(5).deconv多项式除运算a=123;c=4.0013.0028.0027.0018.00d=deconv(c,a)d=4.005.006.00d,r=deconv(c,a)余数余数c除以除以a后的商后的商
5、conv(a,d)%因余数为零，可通过a,d的乘积返回原多项式ans=413282718Matlab的高级计算(6).多项式微分matlab提供了polyder函数多项式的微分。命令格式：polyder(p):求p的微分polyder(a,b):求多项式a,b乘积的微分p,q=polyder(a,b):求多项式a,b商的微分例：a=12345;poly2str(a,x)ans=x4+2x3+3x2+4x+5b=polyder(a)b=4664poly2str(b,x)ans=4x3+6x2+6x+4Matlab的高级计算 k = conv(p,q)k,r = deconv(p,q) k = p
6、olyder(p) k = polyder(p,q)k,d = polyder(p,q) y = polyval(p,x) x = roots(p)多项式运算小结多项式运算中，多项式运算中，使用的是多项式使用的是多项式 系数向量系数向量，不涉及符号计算！不涉及符号计算！Matlab的高级计算4.数值微积分数值微积分4.1 数值积分数值积分 数值积分基本原理数值积分基本原理 求解定积分的数值方法多种多样，如简单的求解定积分的数值方法多种多样，如简单的梯梯形法形法、欧拉法、辛普生欧拉法、辛普生(Simpson) 法法、牛顿牛顿柯特斯柯特斯(Newton-Cotes)法等法等 基本思想都是将整个积分
7、区间基本思想都是将整个积分区间a,b分成分成n个子个子区间区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。这样求定积分问题就分解为求和问题。Matlab的高级计算被积函数由一个表格定义被积函数由一个表格定义 在在MATLAB中，对由表格形式定义的函数关系中，对由表格形式定义的函数关系的求定积分问题用的求定积分问题用trapz(X,Y)函数。其中向量函数。其中向量X,Y定义函数关系定义函数关系Y=f(X)。例例 用用trapz函数计算定积分。函数计算定积分。命令如下：命令如下：X=1:0.01:2.5;Y=exp(-X); %生成函数关系数
8、据向量生成函数关系数据向量trapz(X,Y)ans = 0.28579682416393Matlab的高级计算例：用数值积分法求在x=0到x=10之间所围面积，并讨论步长和积分方法对精度的影响建模：将x的被积区间分为n段，各段长度为xi（i=1，2，n+1）,欧拉法数值积分为梯形法为Matlab的高级计算MATLAB程序fordx=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;%求取样点上的y1n=length(x1);s=sum(y1(1:n
9、-1)*dx;%用欧拉法求积分q=trapz(y1)*dx;%用梯形法求积分stairs(x1,y1),plot(x1,y1)%画出欧拉法及梯形法的积分区域dx,s,q,pause(1),holdoff%给出不同步长下的结果endMatlab的高级计算运行结果：步长欧拉法解梯形法解291081018658150.5841.25816.250.1821.65816.65解析法得到的精确解为2450/3=816.6667右图显示曲线的积分面积，在曲线斜率为负的情况下，欧拉法解偏大，梯形法偏小，精确解位于两者之间步长相同时，梯形法的精度比欧拉法要高Matlab的高级计算变步长辛普生法变步长辛普生法
10、基于变步长辛普生法，基于变步长辛普生法，MATLAB给出了给出了quad函数来求定积分。该函数来求定积分。该函数的调用格式为：函数的调用格式为： I,n=quad(fname,a,b,tol,trace) 其中其中fname是被积函数名。是被积函数名。a和和b分别是定积分的下限和上限。分别是定积分的下限和上限。 tol用用来控制积分精度，缺省时取来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过控制是否展现积分过程，若取非程，若取非0则展现积分过程，取则展现积分过程，取0则不展现，缺省时取则不展现，缺省时取trace=0。 需要事先建立函数文件需要事先建立函数文件fnam
11、e上例中，建立函数文件ex5f.mfunctiony=ex5f(x)y=-x.*x+115;应用quad函数s=quad（ex5f,0,10）运行结果：s=816.6667Matlab的高级计算牛顿柯特斯法牛顿柯特斯法 基于牛顿柯特斯法，基于牛顿柯特斯法，MATLAB给出了给出了quad8函数来求定积分。该函数的调用格函数来求定积分。该函数的调用格式为：式为： I,n=quad8(fname,a,b,tol,trace)其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。 该函数可以更精确地该函数可以更精确地求出定积分的值，且一般情况下函数调
12、求出定积分的值，且一般情况下函数调用的步数明显小于用的步数明显小于quad函数，从而保证函数，从而保证能以更高的效率求出所需的定积分值。能以更高的效率求出所需的定积分值。Matlab的高级计算例例： 分别用分别用quad函数和函数和quad8函数求定积分的函数求定积分的 近似值，并在相同的积分精度下，比较函数近似值，并在相同的积分精度下，比较函数 的调用次数。的调用次数。调用函数调用函数quad求定积分：求定积分：format long;fx=inline(exp(-x); %定义函数定义函数I,n=quad(fx,1,2.5,1e-10)结果：结果：I = 0.28579444254766n
13、 = 65Matlab的高级计算调用函数调用函数quad8求定积分：求定积分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)运行结果运行结果I = 0.28579444254754n = 33Matlab的高级计算二重定积分的数值求解二重定积分的数值求解 使用使用MATLAB提供的提供的dblquad函数就可函数就可以直接求出上述二重定积分的数值解。以直接求出上述二重定积分的数值解。 该该函数的调用格式为：函数的调用格式为： I=dblquad(f,a,b,c,d,tol,trace) 该函数求该函数求f(x,y)在在a,bc,d
14、区域上的二重定区域上的二重定积分。积分。 参数参数tol，trace的用法与函数的用法与函数quad完全相同。完全相同。Matlab的高级计算例例 计算二重定积分计算二重定积分(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.57449318974494Matlab的高级计算4.2 数值微分数值微分 在在MATLAB中，没有直接提供求数值导中，没有直接提供求数值导数的函数，只有计算向前差分的函数
15、数的函数，只有计算向前差分的函数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，按行计算差分。，按行计算差分。Matlab的高级计算例例 生成以向量生成以向量V=1,2,3,4,5,6为基础的范得蒙矩为基础的
16、范得蒙矩阵，按列进行差分运算。阵，按列进行差分运算。命令如下：命令如下：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 0Matlab的高级计算5.数据插值Matlab的高级计算插值
17、方法一维插值的定义已知n个节点，求任意点处的函数值。分段线性插值 linear多项式插值 cubic样条插值 spline最近邻插值方法 nearesty=interp1(x0,y0,x,method)Matlab的高级计算举例举例 多项式插值原理多项式插值原理已已 知知 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)Matlab的高级计算一维插值一维插值: yi = interp1(x, y, x
18、i, 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,yi)Matlab的高级计算Matlab的高级计算%已知数据t=1900:10:1990;p=75.99591.972105.711123.203131.669.150.697179.323203.212226.505249.633;%使用不同
19、的方法进行插值运算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);holdon;plot(x,yi_linear,g,LineWidth,1.5);holdonplot(x,yi_spline,y,LineWidth,1.5);gridon;title(LinearVsSpline);subplot(
20、2,1,2);plot(t,p,ko);holdon;plot(x,yi_cubic,m,LineWidth,1.5);holdonplot(x,yi_nearest,k,LineWidth,1);gridon;title(CubicVsnearest);对某城市人口数量统计统计插值Matlab的高级计算Matlab的高级计算二维插值二维插值高维插值的一种，主要应用于图像处理和数据的可视化，对两个变量的函数z=f(x,y)进行插值zi=interp2(x, y, z, xi, yi, method) method 用法类似一维插值，只是参数用法类似一维插值，只是参数linear表示的是双线性插
21、值方法表示的是双线性插值方法Matlab的高级计算例：以二元函数为例，首先获取基础数据，而后使用二维插值得到更细致的数据，完成绘图%生成基础测量数据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向量数据变为矩阵数据，生成网格点
22、Matlab的高级计算其他插值方法：除了MATLAB内置的插值函数以外，还可以根据工程实践的需要，编写适用于具体领域的插值函数文件来加以调用，如：牛顿插值方法拉格朗日（Lagrange）插值Chebyshev多项式插值方法Matlab的高级计算6.多项式拟和在科学与工程领域，曲线拟合的主要功能是寻求平滑的曲线来最好的表现带有噪声的测量数据，从测量数据中找到变量函数的变化趋势，最后得到曲线拟合的函数表达式Matlab的高级计算已知数据表已知数据表 x x1 x2 xmf(x) y1 y2 ym求拟合函数求拟合函数: (x) = a0 + a1x + + anxn使得使得离散数据的多项式拟合离散数
23、据的多项式拟合达到最小达到最小Matlab的高级计算设i=1,2,N表示按拟合直线求得的近似值，一般地说，它不同于实测值，两者之差称残差。显然，残差的大小是衡量拟合好坏的重要标志，具体地说，可以采用下列三种准则：使残差的最大绝对值为最小：使残差的绝对值之和最小：使残差的平方和为最小：直线拟合直线拟合Matlab的高级计算对于给顶的数据点直线拟合问题可用数学语言描述如下：,求作一次式,使总误差最小。Matlab的高级计算有时候所给出数据点用直线拟合不合适，这时可考虑用多项式拟合，用数学语言描述如下：对于给定的一组数据，寻求作次多项式（）使总误差为最小。曲线拟合曲线拟合Matlab的高级计算多项式
24、拟合命令多项式拟合命令P=polyfit(X,Y,n)求出求出(最小二乘意义下最小二乘意义下)n次拟合多项式次拟合多项式P(x)=a0xn + a1xn-1 + + an-1x + an计算结果为系数计算结果为系数: P= a0, a1, , an-1, an 多项式求值命令多项式求值命令:y1=polyval(P, x)其中其中,P是是n次多项式的系数次多项式的系数,x是自变量的值是自变量的值,y1是是多项式在多项式在x处的值处的值Matlab的高级计算例：如何预报人口的增长例：如何预报人口的增长人口的增长是当前世界上引起普遍关注的问题，并且我们会发现在不同的刊物预报同一时间的人口数字不相同
25、，这显然是由于用了不同的人口模型计算的结果。例如：1949年1994年我国人口数据资料如下：年份xi1949195419591964196919741979198419891994人口数yi5.46.06.77.08.19.19.810.311.311.8建模分析我国人口增长的规律,预报1999年我国人口数。Matlab的高级计算介绍两个简单模型介绍两个简单模型模型一：假设：人口随时间线性地增加模型一：假设：人口随时间线性地增加拟合的精度: Q = ei 2 = (yi - a b xi)2,误差平方和。模型：y = 1.93 + 0.146 x参数估计 观测值的模型： yi = a + b
26、xi + ei ，i = 1,n模型：y = a + b x可以算出：a = 1.93， b = 0.146Matlab的高级计算模型二：模型二： 指数增长模型指数增长模型（用简单的线性最小二乘法）（用简单的线性最小二乘法）用用Matlab软件计算得：软件计算得：a=2.33,b=0.0179即即:Matlab的高级计算程序如下：程序如下：x=1949195419591964196919741979198419891994;y=5.46.06.77.08.19.19.810.311.311.8;a=polyfit(x,y,1);x1=1949:10:1994;y1=a(2)+a(1)*x1;b
27、=polyfit(x,log(y),1);y2=exp(b(2)*exp(b(1)*x1);plot(x,y,*)holdonplot(x1,y1,-r)holdonplot(x1,y2,-k)legend(原曲线,模型一曲线,模型二曲线)Matlab的高级计算程序执行后得到下面图形程序执行后得到下面图形Matlab的高级计算结论的比较如下表：结论的比较如下表：年 份 xi1949195419591964196919741979198419891994人口数 yi5.466.778.19.19.810.311.311.8模型一值5.245.976.77.438.168.99.6210.3611
28、.0911.82误 差0.160.030-0.43-0.060.20.18-0.060.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.51Matlab的高级计算结果分析：结果分析：模型模型 I 2005年年13.43亿，亿， 2010年年14.16亿亿1.误差误差Q1 = 0.2915 m时，此方程成为“超定”方程当n A=1 2 1; 1 0 1;
29、1 3 0; b=2;3;8; x=linsolve(A,b)b是列向量！是列向量！Matlab的高级计算非线性方程的根q Matlab 非线性方程的数值求解非线性方程的数值求解fzero(f,x0)：求方程求方程 f=0 在在 x0 附近的根。附近的根。l 方程可能有多个根，但方程可能有多个根，但 fzero 只给出距离只给出距离 x0 最近的一个最近的一个l fzero 先找出一个包含先找出一个包含 x0 的区间，使得的区间，使得 f 在这个区间在这个区间两个端点上的函数值异号，然后再在这个区间内寻找方程两个端点上的函数值异号，然后再在这个区间内寻找方程 f=0 的根；如果找不到这样的区间
30、，则返回的根；如果找不到这样的区间，则返回 NaN。l x0 是一个标量，不能缺省是一个标量，不能缺省Matlab的高级计算非线性方程的根q 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 不是方程！也
31、不能使用符号表达式！不是方程！也不能使用符号表达式！Matlab的高级计算例： 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 求零点时可以先通过作图确定零点的大致范围求零点时可以先通过作图确定零点的大致范围Matlab的高级计算8.傅立叶分析傅立叶变换是分析周期或非周期信号的频率特性的数学工具，在信号处理领域有着广泛应用Matlab的高级
32、计算例：分析周期信号其中w1=50HZ,w2=120HZt=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)Matlab的高级计算%进行傅立叶变换Y=fft(y,512);%512点的fftPyy=Y.*conj(Y)/512;%conj复共轭函数f=1000*(0:256)/512;%绘图figure,plot(f,Pyy(1:257)Matlab的高级计算N=find(Pyy20);%找出频率较弱的信号Pyy(N)=0;%去除figure,plot(f,Pyy(1:257)Matlab的高级计算
