第四章数值计算ppt课件

1、第四章第四章 数值计算数值计算数值计算内容数值计算内容:数值计算包括数值计算包括:多项式运算多项式运算, 线性方程组求解线性方程组求解, 矩阵特征值问题矩阵特征值问题的解的解, 卷积卷积, 数据分析数据分析, 泛函指令的运用泛函指令的运用, 信号处置和系统分析信号处置和系统分析.4.1 多项式运算多项式运算多项式运算是数学中最根本的运算之一多项式运算是数学中最根本的运算之一.多项式普通可多项式普通可以表示为以表示为: f(x)=a0 xn+a1xn-1+a2xn-2+an-1x+an对于这种表示方式，很容易用一个行向量表示，即：对于这种表示方式，很容易用一个行向量表示，即： T=a0,a1,a

2、2,an-1,an在在MATLAB中中, 多项式正是用这样一个行向量表示的多项式正是用这样一个行向量表示的,系数是按降序陈列的系数是按降序陈列的幂幂4.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+4r=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 = sym4.1.2 多项式的运算多项式的运算多项式的运算主要包括多项式的四那么运算多项式的运算主要包括多项式的四那么运算,导数运导数运算算,估值运算估

4、值运算,求根运算以及多项式的拟合等求根运算以及多项式的拟合等1 多项式的四那么运算多项式的四那么运算多项式四那么运算主要是多项式的加多项式四那么运算主要是多项式的加,减减,乘乘,除除.其中其中,多多项式的加减运算要求两个相加项式的加减运算要求两个相加,减的多项式的大小必需减的多项式的大小必需相等相等,此时加此时加,减才有效减才有效. 当两个相加当两个相加,减的多项式阶次不减的多项式阶次不同时必需经过补同时必需经过补0使其一样使其一样.加加/减减-+/-乘/除-conv, deconvT=deconv(T1, T3)T, r=deconv(T1,T3)商多项式商多项式余式余式T3为分母T1=2,

5、 5, 0, 4, 1, 4;T2=0, 0, 5, 1, 3, 2;T3=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)=

6、5x3+x2+3x+2例例4.1.1-4 多项式求值多项式求值,求上式求上式f1(x)在在x=0.5处的函数值处的函数值T1=2,5,0,4,1,4;x=0.5;y=polyval(T1,x)y = 5.87503 多项式的导数运算多项式的导数运算-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,

7、4; root=roots(T1);root= -2.7709 0.5611 + 0.7840i 0.5611 - 0.7840i -0.4257 + 0.7716i -0.4257 - 0.7716i10*x4+20*x3+8*x+14 拟合和插值拟合和插值-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,

8、s = polyfit(x0, y0, n);x0,y0-给定数据对给定数据对n-拟合出的多项式次数拟合出的多项式次数p-多项式向量多项式向量s-偏向信息偏向信息yi=interp1(x0, y0, xi, cubic);xi, yi-得到的新的插值点得到的新的插值点cubic-插值方式插值方式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;n=3; % 设定拟合次数为设定拟合次数为3p,s=polyfit(x0,y0,n); % 得到拟合多项式向量和相关偏向信得到拟合多项式向量和相关偏向信息息 x

9、x=0:0.01:1; yy=polyval(p,xx); % 按拟合曲线计算采样值按拟合曲线计算采样值n1=6; % 设定拟合次数为设定拟合次数为6p1,s1=polyfit(x0,y0,n1);yy1=polyval(p1,xx);plot(xx,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); y=51.48x3 -77.74x2+35.06x-0.20y=348.21x6-1060.23x

10、5+1297.68x4-758.90 x3+ 181.00 x2 + 1.00 x+0.48拟合只能在给定数据所限定的区间内运用拟合只能在给定数据所限定的区间内运用, 不要恣意向外拓展不要恣意向外拓展例例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.45,5.35,9.22;xi=0:0.02:1;yi=interp1( x0, y0, xi, cubic ); subplot(1,2,1);plot(xi, yi

11、, -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);插值插值拟合拟合插值方式插值方式:linear-线性插值线性插值nearest-最近插值最近插值cubic-三次多项式插值三次多项式插值spline-样条插值样条插值pchip-分段三次分段三次Hermite插值插值v5c

12、ubic-MATLAB 5.0版的三次多项式插值版的三次多项式插值.拟合和插值的区别拟合和插值的区别:曲线拟合是研讨如何寻觅曲线拟合是研讨如何寻觅“平滑平滑曲线最好地表现带噪声的曲线最好地表现带噪声的“丈量数据丈量数据,但并不要求拟但并不要求拟合曲线经过这些合曲线经过这些“测试数据点测试数据点.插值是在认定所给数插值是在认定所给数据完全正确的情况下据完全正确的情况下,研讨如何平滑地估算出研讨如何平滑地估算出“基准数基准数据之间其他点的函数值据之间其他点的函数值,因此插值一定经过因此插值一定经过“基准数基准数据据通俗地讲通俗地讲, 拟合就是由知点得到一条曲线拟合就是由知点得到一条曲线, 使该曲线

13、最使该曲线最能反映点所代表的规律能反映点所代表的规律.比如做欧姆定理的实验的时候比如做欧姆定理的实验的时候,由于实验中存在误差由于实验中存在误差,最后拟合得到的曲线是一条直线最后拟合得到的曲线是一条直线,而且一定只需部分点落在拟合的直线上而且一定只需部分点落在拟合的直线上,但此时该直线但此时该直线和测试点的方差最小和测试点的方差最小.由拟合直线的斜率就可以知道电由拟合直线的斜率就可以知道电阻的阻值阻的阻值.拟合是探测事物变化规律的方法拟合是探测事物变化规律的方法.插值就是根据函数上某些知点插值就是根据函数上某些知点(或实验数据或实验数据),按一定规按一定规律律(插值方法插值方法)寻求未知的点寻

14、求未知的点,比如知一个常用对数比如知一个常用对数y=log(x)表表,是按照是按照x=0.1:0.1:10制表的制表的,假设按假设按知数据求知数据求y=log(2.897)就可以用插值得到就可以用插值得到.表制得表制得越密越密,插值越准确插值越准确例例4.1.1-9 知经过实验得到电压和电流的数值如下：知经过实验得到电压和电流的数值如下：i0=0:0.5:10;v0=0,98,202,306,395,506,592,708,789,890,1009,1108,1186, 1310,1,1510,1601,1716,1782,1920,2019;按一次按一次,二次二次,三次多项式对数值关系进展拟

15、合三次多项式对数值关系进展拟合阐明：阐明：为保证较好的拟合效果，多项式阶数要获得适当，过为保证较好的拟合效果，多项式阶数要获得适当，过低，残差较大，过高，拟合模型将包含噪声影响，通低，残差较大，过高，拟合模型将包含噪声影响，通常要保证拟合阶数小于数据对的数目。常要保证拟合阶数小于数据对的数目。i0=0:0.5:10;v0=0,98,202,306,395,506,592,708,789,890,1009,1108, . 1186,1310,1,1510,1612,1716,1782,1920,2019;n1=1; p1,s=polyfit(i0,v0,n1); n2=2; p2,s=polyf

16、it(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(二次拟合二次拟合);subplot(1,3,3),plot(ii,v3,-b,i0,v0,.r,MarkerSize,8),

17、 xlabel(三次拟合三次拟合);p1=201.00, -2.89P2=0.31, 197.90, 2.02P3=0.03, -0.10, 199.49, 0.85p1=201.00,-2.89P2=0.31, 197.90, 2.02P3=0.03, -0.10, 199.49, 0.85对于含对于含n个未知数的个未知数的n个方程构成的方程组个方程构成的方程组Ax=b, 在线性代数教在线性代数教科书中科书中,最常引见的解法有最常引见的解法有:Cramer法法;逆阵法逆阵法, 即即x=A-1b;高斯消高斯消元法元法; LU法法 对于维数不高对于维数不高,条件数不大的矩阵条件数不大的矩阵,以上

18、四种所得结果不会有明显以上四种所得结果不会有明显的差别的差别.但前三种解法的更多的意义在实际上但前三种解法的更多的意义在实际上,而不在实践的数值而不在实践的数值计算上计算上.在在MATLAB中中,方程采用方程采用LU法求解法求解,并且出于算法稳定性并且出于算法稳定性的思索的思索,行列式和逆的计算也都是在行列式和逆的计算也都是在LU分解的根底上进展的分解的根底上进展的MATLAB中的四条指令中的四条指令:L, U, P=lu(A)-矩阵的矩阵的LU分解分解,使使LU=PA. 多种格式多种格式.det(A)-求矩阵求矩阵A的行列式的行列式, 它由它由detA=detUuii 算得算得inv(A)-

19、求矩阵求矩阵A的逆的逆, 由由A-1=U-1L-1P算得算得x=Ab-除法指令求方程的解除法指令求方程的解(对恰定方程对恰定方程, 采用采用LU法执行法执行)4.2.1 方程组求解方程组求解4.2 线性方程组的解线性方程组的解对于方程组对于方程组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+7x3=1-2x1+4x2+5x3=-7A=2, 2, 3; 4, 7, 7; -2, 4

20、, 5, b=3; 1; -7x=Ab, y=b/A求解方程时求解方程时,尽量不要运用指令尽量不要运用指令inv(A)*b进展进展,它不仅计算速度它不仅计算速度没有除法快没有除法快,而且计算精度也没除法高而且计算精度也没除法高.此外此外,除法还适用于除法还适用于“超超定和定和“欠定方程欠定方程x = 9.00 -36.00 31.00y = 9.00 -36.00 31.00方程组的解的三种情况方程组的解的三种情况:对于方程对于方程Ax=b, AAx=b, A为为AmAmn n矩阵矩阵, ,有三种情况有三种情况: : 当当m=nm=n时时, ,此方程成为此方程成为 恰定恰定 方程方程, ,求解

21、准确解求解准确解 当当mnmn时时, ,此方程成为此方程成为“超定方程超定方程, ,寻求最小二乘解寻求最小二乘解 ( (直线拟直线拟合合) ) 当当mnmn时此时无解时此时无解,但适用中可以求最小二乘解即但适用中可以求最小二乘解即:方程解方程解 (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;

22、2;3;x1=inv(A*A)*A*b;x2=Abx1 = 1.00 0.00 x2 = 1.00 0 3) 3) 欠定方程欠定方程当方程数少于未知量个数时当方程数少于未知量个数时, ,即不定情况即不定情况, ,有无穷多个解存在有无穷多个解存在. .matlab可求出两个解：可求出两个解：a.用除法求的解用除法求的解x是具有最多零元素的解是具有最多零元素的解b.具有最小长度或范数的解具有最小长度或范数的解,这个解是基于伪逆这个解是基于伪逆pinv求得的求得的.例例4.2.1-4求以下欠定方程组的解求以下欠定方程组的解x1+2x2+3x3=12x1+3x2+4x3=2A=1,2,3;2,3,4;

23、b=1;2;x1=Ab;x2=pinv(A)*bx1 = 1 0 0 x2 = 0.8333 0.3333 -0.16674.2.2 奇特值分解和矩阵构造奇特值分解和矩阵构造相应的相应的MATLAB指令如下指令如下:U,S,V=svd(A)-矩阵的奇特值分解三对组阵，使矩阵的奇特值分解三对组阵，使A=USVTnorm(A,flag)-计算矩阵计算矩阵A的范数，范数类型由的范数，范数类型由flag指定指定.cond(A)-计算矩阵计算矩阵A的条件数的条件数.pinv(A)-求矩阵的广义逆求矩阵的广义逆flag-1, 2, inf, fro4.2.3 线性二乘问题的解线性二乘问题的解对于超定方程对

24、于超定方程,求其最小二乘解有三种方法：求其最小二乘解有三种方法：正那么方程法得解正那么方程法得解x=(A A)-1A b, 数值精度差数值精度差广义逆法得解广义逆法得解x=A+b, 所得解可靠所得解可靠.用矩阵除法得解用矩阵除法得解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)

25、*y, xxx=Ayx = 3.4811 5.1595 0.9534 -0.0466xx = 3.4759 5.1948 0.7121 -0.1101xxx = 3.4605 5.2987 0 -0.2974由于阵缺一个列秩由于阵缺一个列秩,除法给出的最小二乘除法给出的最小二乘解就只需三个非零元解就只需三个非零元素素.它是只需三个非零它是只需三个非零元素的一切最小二乘元素的一切最小二乘解中范数最小的解中范数最小的Gallery-MATLAB 设置设置的特殊矩阵的特殊矩阵. (2) 计算三个解的范数计算三个解的范数nx=norm(x), nxx=norm(xx), nxxx=norm(xxx)n

26、x = 6.2968nxx = 6.2918nxxx = 6.3356用广义逆所得的最小二乘范数最小(3) 比较三种解法的方程误差比较三种解法的方程误差e=norm(y-A*x), ee=norm(y-A*xx), eee=norm(y-A*xxx)e = 0.6986ee = 0.0474eee = 0.0474误差的平方根最大，精度最差误差的平方根最大，精度最差4.2.4 特征值分解和矩阵函数特征值分解和矩阵函数涉及矩阵特征值分解的常用的涉及矩阵特征值分解的常用的MATLAB指令如下指令如下:d = eig(A)-计算矩阵计算矩阵A的特征值的特征值,以向量方式存放以向量方式存放V,D=ei

27、g(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, D=eig(A)和和c=condeig(A)两两条指令条指令例例4.2.4-1 验

28、证验证det(A)=iA=2,-1,-1; 3,4,-2; 3,-2,4det(A)d=eig(A)ans = 60d = 2.0000 + 2.4495i 2.0000 - 2.4495i 6.0000 d(1)*d(2)*d(3)ans = 60.0000行列式行列式特征值特征值% 矩阵行列式等于其一切特征值的积矩阵行列式等于其一切特征值的积4.3 数据分析函数数据分析函数(1) 假设输入宗量假设输入宗量x是向量是向量,那么不论是行向量还是列向量那么不论是行向量还是列向量,运算是运算是对整个向量进展的对整个向量进展的(2) 假设输入宗量假设输入宗量x是二维数组，那么指令运算是按列进展的是二

29、维数组，那么指令运算是按列进展的MATLAB在进展数据分析时的商定在进展数据分析时的商定:4.3.1 随机数发生器和统计分析指令随机数发生器和统计分析指令函数函数含义含义rand(m, n)产生产生mxn维的维的0,1区间均匀分布随机数组区间均匀分布随机数组randn(m,n)产生产生mxn维的均值为维的均值为0标准差为标准差为1的正态分布随的正态分布随机数组机数组min(x)对矩阵各列分别求最小值对矩阵各列分别求最小值max(x)对矩阵各列分别求最大值对矩阵各列分别求最大值median(x)对矩阵各列分别求中位数对矩阵各列分别求中位数mean(x)对矩阵各列分别求均值对矩阵各列分别求均值st

30、d(x)对矩阵各列分别求标准差对矩阵各列分别求标准差var(x)对矩阵各列分别求方差对矩阵各列分别求方差C=conv(x)给出矩阵各列的协方差阵给出矩阵各列的协方差阵P=corrcoef(x) 给出矩阵各列间的相关系数给出矩阵各列间的相关系数例例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, c

31、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 94.3.2 差分和累计指令差分和累计指令 函数函数含义含义del2(U,hx,hy)五点五点Laplaciandiff(X, m, n)沿第沿第n维求维求m次差分次差分 (二维：二维：n=1; 列列,n=2,行行)DZx,DZy=grad

32、ient(z,dx,dy)对对Z求求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元素按模增大排列元素按模增大排列diff, trapz, cumtrapz指令的演示指令

33、的演示a=1,2,3,4,5; b=a, a+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.50 2.50 3.50 4.504.00 6.00 8.00 10.007.50 10.50 13.50 16.5012.00 16.00 20.00 24.001 2 3 42 3 4 53 4 5 64 5 6 75 6 7 8例例4.3.2-1求求1+2+3+100以及以及50! x=1:1:100;s

34、um_x=sum(x);sum_x = 5050a=1:1:50;prod_a=prod(a)prod_a = 3.0414e+064cumsum, cumprod的用法的用法上两个函数分别是求向量的累计和和累计乘积上两个函数分别是求向量的累计和和累计乘积.假设假设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*n例例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*tra

35、pz(y)s=dt*cumsum(y);s3=s(end)s=dt*cumtrapz(y);s4=s(end)matlab还有更精良的积分指令还有更精良的积分指令: quad, quda8s1 = 8.0060s2 = 8.000001s3 = 8.0060s4 = 8.0000014.4 MATLAB泛函指令泛函指令在在MATLAB中中,凡以函数为输入宗量的指令凡以函数为输入宗量的指令,都被统称为泛函指令都被统称为泛函指令. 最常见的有最常见的有:z=fzero(fun,x0)-求一元函数零点指令的最简单格式求一元函数零点指令的最简单格式x=fsolve(fun,x0)-解非线性方程组的最简

36、单格式解非线性方程组的最简单格式x=fminbnd(fun,x1,x2)-求函数在区间求函数在区间(x1,x2)中极小值的指令中极小值的指令 最简格式最简格式x=fminsearch(fun,x0)-单纯形法求多元函数极值点指令的最单纯形法求多元函数极值点指令的最简格式简格式.x=fminunc(fun,x0)-拟牛顿法求多元函数极值点指令的最简拟牛顿法求多元函数极值点指令的最简格式格式a=lsqnonlin(fun,a0)-解非线性最小二乘问题指令的最简格式解非线性最小二乘问题指令的最简格式.q=quad(fun,a,b)-采用递推自顺应采用递推自顺应simpson法计算积分法计算积分q=q

37、uadl(fun,a,b)-采用递推自顺应采用递推自顺应Lobatto法求数值积分法求数值积分SS=dblquad(fun,inmin,inmax,outmin,outmax)-二重闭环数二重闭环数值积分值积分t,YY=ode45(fun,tspan,Y0)-采用采用4,5阶阶Runge-Kutta方程解方程解算算ODE初值问题初值问题指令中被处置的函数指令中被处置的函数fun,可以取可以取:字符串表达式字符串表达式,内联函数内联函数, M函数函数文件的函数句柄文件的函数句柄4.4.1 求函数零点求函数零点对于恣意函数对于恣意函数f(x)=0,零点情况复杂零点情况复杂,没有一个通用解法没有一个

38、通用解法.普普通来说通来说,零点的数值计算过程是零点的数值计算过程是:先猜测一个初始零点或先猜测一个初始零点或该零点所在的区间该零点所在的区间,然后经过计算然后经过计算,使猜测值不断准确化使猜测值不断准确化,或使猜测区间不断收缩或使猜测区间不断收缩,直到到达预定的精度直到到达预定的精度,终止计算终止计算.求函数的零点的步骤为求函数的零点的步骤为:(1) 利用利用matlab作图指令获取初步近似值作图指令获取初步近似值详细做法详细做法: 先确定一个零点能够存在的自变量区间先确定一个零点能够存在的自变量区间,然后利用然后利用plot指令画出指令画出f(x)在该区间中的图形在该区间中的图形,察看察看

39、f(x)与横轴的交点坐标与横轴的交点坐标,也可以也可以用用zoom对交点处部分放大再读数对交点处部分放大再读数.借助借助ginput指令获得更准确的指令获得更准确的交点坐标值交点坐标值(2) 利用利用fzero指令求准确解指令求准确解z=fzero(fun,x0)-求一元函数零点指令的最简格式求一元函数零点指令的最简格式z, f_z, exitflag=fzero(fun,x0,options,p1,p2,)-最完好格式最完好格式 阐明：阐明： 由于由于fzero是根据函数能否越横轴来决议零点的是根据函数能否越横轴来决议零点的,因此本程序无法确定函数因此本程序无法确定函数曲线仅触及横轴和不穿越

40、的零点曲线仅触及横轴和不穿越的零点 第二个输入量第二个输入量x0表示零点的初始猜测表示零点的初始猜测.他可以是标量或二元向量他可以是标量或二元向量.当当x0取标取标量时量时,该指令将在它两侧寻觅一个与之最接近的零点该指令将在它两侧寻觅一个与之最接近的零点;当取二元向量当取二元向量a,b时时,该该指令将在区间指令将在区间a,b内寻觅一个零点内寻觅一个零点. 输入宗量输入宗量options是优化迭代采用的参数是优化迭代采用的参数.P1,P2是向函数是向函数fun传送的参数传送的参数. z是所求零点的自变量值是所求零点的自变量值. f_z是函数值是函数值,exitflag阐明程序终止计算的条件阐明程

41、序终止计算的条件.假设假设exitflag0, 阐明找到零点后退出阐明找到零点后退出例例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); % 便于查看零便于查看零点点xlabel(x), ylabel(y);title(曲线方程式为：曲线方程式为：y=2x-10*x);(2) 获取初步近似解获取初步近似解zoom on;tt,yy=ginput(2);zoom off;fun

42、c=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);4.4.2 求函数极值点求函数极值点许多科学研讨和工程计算问题都可以归结为一个极值问题许多科学研讨和工程计算问题都可以归结为一个极值问题,如能量最小，时如能量最小，时间最短间最短,最正确拟合最正确拟合,最短途径等最短途径等.

43、又由于又由于f(x)的极小值问题等价于的极小值问题等价于-f(x)的极大的极大值问题值问题,所以所以matlab只需处置极小值的指令只需处置极小值的指令确切地说确切地说,这里讨论的只是这里讨论的只是“局域极值问题局域极值问题.“全域最小问题全域最小问题要复杂得多要复杂得多.至今没有一个系统性的方法可求解普通的至今没有一个系统性的方法可求解普通的“全域最全域最小问题小问题.对于一元二元函数对于一元二元函数,作图察看对确定全域最小值有很作图察看对确定全域最小值有很好的运用价值好的运用价值.但更多元的函数但更多元的函数,就很难运用作图法就很难运用作图法.matlab求函数极值的三条指令求函数极值的三

44、条指令x, fval, exitflag, output=fminbnd(fun,x1,x2,opt,p1,p2,)-求一元函数在区间求一元函数在区间(x1,x2)中极小值的指令的最完好格式中极小值的指令的最完好格式x, fval, exitflag, output=fminsearch(fun,x0,opt,p1,p2,)-单纯形法求多元函数极值点指令的最完好格式单纯形法求多元函数极值点指令的最完好格式x,fval,exitflag,output,grad,hessian=fminunc(fun,x0,opt,p1,p2,)-拟牛顿法求多元函数极值点指令的最完好格式拟牛顿法求多元函数极值点指

45、令的最完好格式例例4.4.2-1 求例求例4.4.1-1 中的函数的最小值中的函数的最小值(1) 函数采用字符串表达式函数采用字符串表达式func=2x-10*x;x1=0,x2=8; % 经过察看经过察看, 发现最小值发现最小值在此区间在此区间Xmin,Ymin=fminbnd(func,x1,x2)Xmin = 3.8507Ymin = -24.0800(2) 函数采用函数采用M文件函数的句柄文件函数的句柄% M文件内容，保管为文件内容，保管为func.mfunction y=f(x)y=2x-10*x; % M文件内容文件内容x1=0; x2=8;Hf=func;Xmin,Ymin=fminbnd(Hf,x1,x2)4.4.3 数值积分数值积分数值积分有闭型和开型算法之分两者的区别在于：能否需求计数值积分有闭型和开型算法之分两者的区别在于：能否需求计算积分区间端点处的函数值算积分区间端点处的函数值matlab中仅提供闭型数值积分指令中仅提供闭型数值积分指令:q=quad(fun, a, b, tol, trace, p1, p2, )-采用递推自采用递推自 顺应顺应Simpson法计算积分法计算积分q=qudal(fun, a, b, tol, trace, p1, p2,)-采用递推自采用递推自顺应顺应Lobatlo法求数值积分