第五章数值计算

1、 具体知识点见导图具体知识点见导图5。5.1 矩阵代数的计算矩阵代数的计算 矩阵代数计算包括求矩阵的各种三角分解、线性方程组的求解、矩阵代数计算包括求矩阵的各种三角分解、线性方程组的求解、矩阵的逆和广义逆、矩阵的模及条件数、矩阵的谱分解和奇异值分矩阵的逆和广义逆、矩阵的模及条件数、矩阵的谱分解和奇异值分解。矩阵代数的计算是解决数学问题的基础性方法，任何的数学模解。矩阵代数的计算是解决数学问题的基础性方法，任何的数学模型不论是统计模型、微分方程模型还是其他的复杂数学模型最后大型不论是统计模型、微分方程模型还是其他的复杂数学模型最后大部分将进行离散化处理，然后归为矩阵代数的计算。部分将进行离散化处

2、理，然后归为矩阵代数的计算。5.1.1 矩阵代数知识点搜索矩阵代数知识点搜索 在参考文献在参考文献Reference中进入中进入MATLAB Function Reference，即进入，即进入MATLAB的函数参考列表，然后双击的函数参考列表，然后双击Mathematics，我们就进入了，我们就进入了MATLAB的核心数学方法的函数的核心数学方法的函数集（或称数学方法命令集）。集（或称数学方法命令集）。 代数计算函数代数计算函数【例【例5.1】MATLAB在数值计算方面的强大功能。在数值计算方面的强大功能。传统求逆法和传统求逆法和MATLAB左除法求线性方程的性能对比。左除法求线性方程的性能

3、对比。rand(state,12);%选定随机种子。选定随机种子。A=rand(100,100)+1.e8; %生成生成100100均匀分布随机矩阵。均匀分布随机矩阵。x=ones(100,1); %令解向量令解向量 x 为全为全1的的100元列向量。元列向量。b=A*x; %为使为使 Ax=b 方程一致，用方程一致，用A和和 x 生成生成 b 向量。向量。cond(A) %求求 A 阵的条件数。阵的条件数。 ans = 1.4426e+012 从矩阵的条件数来看，相当大。说明矩阵从矩阵的条件数来看，相当大。说明矩阵A是严重病态矩阵，求以是严重病态矩阵，求以它为系数矩阵的线性方程的根将是不精确

4、的，或者说是不稳定的。它为系数矩阵的线性方程的根将是不精确的，或者说是不稳定的。% “求逆求逆”法解恰定方程的误差、残差、运算次数和所用时间法解恰定方程的误差、残差、运算次数和所用时间flops(0);tic %计数器置计数器置0，启动计时器，启动计时器Stopwatch Timerxi=inv(A)*b; % xi 是用是用“求逆求逆”法解恰定方程所得的解。法解恰定方程所得的解。ti=toc %关闭计时器，并显示解方程所用的时间。关闭计时器，并显示解方程所用的时间。ci=flops%“求逆求逆”法解方程所用的运算次数法解方程所用的运算次数eri=norm(x-xi) %解向量解向量 xi 与

5、真解向量与真解向量 x 的范的范-2误差。误差。rei=norm(A*xi-b)/norm(b) %方程的范方程的范-2相对残差相对残差 ti = 0.9300 ； 用时用时0.93秒秒ci = 2070322 ； 乘除次数乘除次数eri = 3.0708e-004 绝对误差绝对误差rei = 6.6280e-007 相对误差相对误差% “左除左除”法解恰定方程的误差、残差、运算次数和所用时间法解恰定方程的误差、残差、运算次数和所用时间flops(0);tic;xd=Ab; % 以上是用以上是用“左除左除”法解方程所得的解。法解方程所得的解。td=toc,cd=flops,erd=norm(x

6、-xd),red=norm(A*xd-b)/norm(b) td = 0.2200 计算时间计算时间cd = 741872 乘除法的次数乘除法的次数erd = 3.2243e-004 绝对误差绝对误差red = 2.0095e-016 相对误差相对误差对比结果：用左除法是逆阵法运算次数的对比结果：用左除法是逆阵法运算次数的4/10，时间是逆阵除法，时间是逆阵除法的的1/4。左除法的相对残差几乎为。左除法的相对残差几乎为 “机器零机器零” ，比逆阵法相对精，比逆阵法相对精度要高的多。度要高的多。5.1.2 矩阵分析矩阵分析Matrix Analysis矩阵分析主要是考察一个系数矩阵矩阵分析主要是

7、考察一个系数矩阵A对线性方程组解的灵敏度，即对线性方程组解的灵敏度，即解的稳定性问题。例如对于最小二乘估计的解解的稳定性问题。例如对于最小二乘估计的解 如果协方差矩阵具有列相关性时，或者说该矩阵接近不满秩如果协方差矩阵具有列相关性时，或者说该矩阵接近不满秩时，我们称该矩阵为病态矩阵，它将严重影响最小二乘估计的稳定时，我们称该矩阵为病态矩阵，它将严重影响最小二乘估计的稳定性。度量其病态性的指标为条件数，矩阵的条件数为该矩阵的最大性。度量其病态性的指标为条件数，矩阵的条件数为该矩阵的最大特征值比最小特征值，用来衡量矩阵计算的精确度，当条件数为特征值比最小特征值，用来衡量矩阵计算的精确度，当条件数为

8、1时是很好的条件，而条件数非常大时则表明矩阵不满秩，计算时将时是很好的条件，而条件数非常大时则表明矩阵不满秩，计算时将产生很大的误差。产生很大的误差。矩阵分析包括矩阵分析包括:条件数的计算、矩阵模的计算、矩阵的行列式的计算和矩阵的标准条件数的计算、矩阵模的计算、矩阵的行列式的计算和矩阵的标准正交基的计算等等。正交基的计算等等。 yXXXb1 【例【例5.1.1】求矩阵】求矩阵A的条件数、行列式、模等的计算的条件数、行列式、模等的计算 7518244518AA=8,-1,-5；-4,4,-2；18,-5,7 % 定义一个矩阵定义一个矩阵c=cond(A) % 计算矩阵计算矩阵A的条件数的条件数c

9、 = 8.2546n = norm(A) % 计算矩阵的模计算矩阵的模n = 21.5249d = det(A) % 计算矩阵的行列式计算矩阵的行列式d = 412B = orth(A) % 求矩阵的一个标准正交基求矩阵的一个标准正交基 BB =B = -0.2987 -0.9498 -0.0933 -0.2987 -0.9498 -0.0933 0.2469 -0.1713 0.9538 0.2469 -0.1713 0.9538 -0.9219 0.2619 0.2856 -0.9219 0.2619 0.2856EYE=B*B % 验证B是标准正交基 EYE =EYE = 1.0000

10、-0.0000 0.0000 1.0000 -0.0000 0.0000 -0.0000 1.0000 0 -0.0000 1.0000 0 0.0000 0 1.0000 0.0000 0 1.00005.1.3 矩阵的分解和三角分解矩阵的分解和三角分解Factorization 从计算的角度来看，处理矩阵问题的一个最有效的方法是，将从计算的角度来看，处理矩阵问题的一个最有效的方法是，将一个一般的矩阵分解为几个简单矩阵的乘积形式。上三角矩阵为例，一个一般的矩阵分解为几个简单矩阵的乘积形式。上三角矩阵为例，由其性质，两个上三角阵的和、积仍为上三角阵，上三角阵的特征由其性质，两个上三角阵的和、积

11、仍为上三角阵，上三角阵的特征值就是其对角线元素。系数矩阵为上三角阵的线性方程组只需用回值就是其对角线元素。系数矩阵为上三角阵的线性方程组只需用回代的方法求解，上三角阵的逆阵仍然是上三角阵。代的方法求解，上三角阵的逆阵仍然是上三角阵。LU分解：分解： 将矩阵将矩阵A分解为分解为A = LU，其中，其中L为单位下三角阵，为单位下三角阵，U为上三角阵。为上三角阵。Cholesky分解：将正定对称阵分解成分解：将正定对称阵分解成 A=TT的形式，其中的形式，其中T为为正交变换且为上三角阵。正交变换且为上三角阵。QR分解：分解： 将矩阵分解为将矩阵分解为A=QR的形式，其中的形式，其中Q为一个正交阵，为

12、一个正交阵，R为一个上三角阵。为一个上三角阵。【例【例5.1.2】对矩阵】对矩阵A进行进行LU分解分解L,U = lu(A) % 对矩阵对矩阵A进行进行LU分解分解A = 8 -1 -5 -4 4 -2 18 -5 -7L = 0.4444 0.4231 1.0000 -0.2222 1.0000 0 1.0000 0 0U = 18.0000 -5.0000 -7.0000 0 2.8889 -3.5556 0 0 -0.3846【例【例5.1.3】对对称正定矩阵进行】对对称正定矩阵进行Cholesky分解分解B = 1 -3 2 -3 10 -5 2 -5 6 T = chol(B) %

13、对对称矩阵进行对对称矩阵进行Cholesky分解分解C = T*T % 对对T进行验证进行验证B = 1 -3 2B = 1 -3 2 -3 10 -5 -3 10 -5 2 -5 6 2 -5 6T = 1 -3 2T = 1 -3 2 0 1 1 0 1 1 0 0 1 0 0 1C = 1 -3 2 C = 1 -3 2 -3 10 -5 -3 10 -5 2 -5 6 2 -5 6 5.1.4 求解线性方程组求解线性方程组Linear Equations【例【例5.1.4】求解以下线性方程组】求解以下线性方程组A=1,-3,2;-3,10,-5;2,-5,6 % 定义系数矩阵定义系数矩

14、阵Ab=27,-78,64 % 定义常数定义常数bx=Ab % 求解线性方程组求解线性方程组b1=A*x最后的解为：最后的解为：x=(1 -4 7) 6478276525103231321xxx5.1.5 矩阵的特征值或奇异值及其与特征向量矩阵的特征值或奇异值及其与特征向量Eigenvalues and Singular Values（1）对称矩阵的谱分解）对称矩阵的谱分解,即求矩阵的特征值和特征向量即求矩阵的特征值和特征向量设矩阵设矩阵A是是NN维方阵，则谱分解定义如下：维方阵，则谱分解定义如下：UUUDUAn1其中其中D为对角矩阵，为对角矩阵，U为正交矩阵。对角矩阵的元素称为矩阵为正交矩阵

15、。对角矩阵的元素称为矩阵A的的特征值，矩阵特征值，矩阵U的列向量为对应于矩阵的列向量为对应于矩阵A的特征值的特征向量。的特征值的特征向量。（2）对一般矩阵的奇异值分解）对一般矩阵的奇异值分解设矩阵设矩阵A为为nm阶矩阵，则矩阵的奇异值分解定义为：阶矩阵，则矩阵的奇异值分解定义为：mmmnnnmnVDUA其中其中 U 和和 V 分别为分别为 n 和和 m 阶正交方阵，阶正交方阵，D 为为 nm 阶对角阵，阶对角阵，其非对角元素均为其非对角元素均为 0，D的对角线元素称为矩阵的对角线元素称为矩阵 A 的奇异值。奇的奇异值。奇异值的分解与矩阵的谱分解方法类似。异值的分解与矩阵的谱分解方法类似。求解矩

16、阵求解矩阵A的特征值和特征向量的语法为：的特征值和特征向量的语法为：V, D = eig(A)其中其中V为矩阵其列向量为矩阵为矩阵其列向量为矩阵A的特征向量，的特征向量，D为对角矩阵对角线为对角矩阵对角线元素为矩阵元素为矩阵A的特征值。的特征值。 【例【例5.1.5】求矩阵的特征值和奇异值。】求矩阵的特征值和奇异值。（1）求上面矩阵）求上面矩阵A的特征值和特征向量。的特征值和特征向量。V,D = eig(A) % 求矩阵求矩阵A的特征值和特征向量的特征值和特征向量V = 0.9654 -0.0091 0.2606 0.2211 0.5581 -0.7998 -0.1381 0.8297 0.5

17、408D = 0.0266 0 0 0 2.6151 0 0 0 14.3583（2）求下列矩阵）求下列矩阵B的奇异值及其对应的向量。的奇异值及其对应的向量。B =9 4；6 8； 2 7 % 矩阵矩阵B为长方形矩阵为长方形矩阵U,S,V = svd(B) % 求矩阵求矩阵B的奇异值和对应的两个向量的奇异值和对应的两个向量B = 9 4 6 8 2 7U = -0.6105 0.7174 0.3355 -0.6646 -0.2336 -0.7098 -0.4308 -0.6563 0.6194S = 14.9359 0 0 5.1883 0 0V = -0.6925 0.7214 -0.721

18、4 -0.69255.2 多项式与插值多项式与插值 在数学建模竞赛中或在科学研究中，使用频率最高的数据加工在数学建模竞赛中或在科学研究中，使用频率最高的数据加工技术是对实验数据的拟合，拟合有很多方法如多项式方法、最技术是对实验数据的拟合，拟合有很多方法如多项式方法、最小二乘法和样条插法。当我们对给定的数据利用上面的方法选小二乘法和样条插法。当我们对给定的数据利用上面的方法选择了拟合方法后，对任意给定的自变量点（插值点）我们可以择了拟合方法后，对任意给定的自变量点（插值点）我们可以计算该点的函数值了。计算该点的函数值了。MATLAB约定多项式如下：约定多项式如下： 1121)(nnnnaxaxa

19、xaxP 我们可以用行向量来表示一个多项式。如，我们可以用行向量来表示一个多项式。如，。我们可以进行多项式的乘、除、因子的提取、多项式的简化、求。我们可以进行多项式的乘、除、因子的提取、多项式的简化、求多项式的根等等常规的多项式的计算和多项式的拟合。多项式的根等等常规的多项式的计算和多项式的拟合。注意：如果某系数为注意：如果某系数为0，必须写上以占一个位置。，必须写上以占一个位置。,121 naaaP5.2.1 多项式的创建方法多项式的创建方法（1）直接用行向量创建，如：）直接用行向量创建，如：P=1，2，3，5，6，7 %这创建了一个这创建了一个5阶多项式。阶多项式。 （2）利用命令：）利用

20、命令：P = ploy（A） % 产生多项式系数向量。产生多项式系数向量。这里这里A为方阵，即多项式为方阵，即多项式 P 为为 A 的特征多项式。的特征多项式。【例【例5.2.1】求三阶魔方方阵】求三阶魔方方阵M 的特征多项式。的特征多项式。M=magic(3)PM=poly(M)； % A的特征多项式的特征多项式PPM=poly2str(PM,x) % 将多项式转成习惯的方式将多项式转成习惯的方式RM=roots(PM) % 求求A或者说特征多项式的根或者说特征多项式的根PPA = x3 - 15 x2 - 24 x + 360RM = 15.0000 -4.8990 4.89907632)

21、(345 xxxxxP5.2.2 计算多项式的操作函数与例题计算多项式的操作函数与例题 多项式的运算包括，多项式相乘、多项式相除、求多项式的公多项式的运算包括，多项式相乘、多项式相除、求多项式的公因式等等运算，表因式等等运算，表5.2.1 列出一些基本多项式运算命令列出一些基本多项式运算命令【例【例5.2.2】由给定的根向量求多项式的系数向量。】由给定的根向量求多项式的系数向量。R=-0.5,-0.3+0.4*i,-0.3-0.4*i; % 根向量根向量P=poly(R) % R的特征多项式的特征多项式PR=real(P) % 求求PR的实部的实部PPR=poly2str(PR,x) 计算结果

22、如下：计算结果如下：P = 1.0000 1.1000 0.5500 0.1250PR = 1.0000 1.1000 0.5500 0.1250PPR = x3 + 1.1 x2 + 0.55 x + 0.125 注意：注意：（1）要形成实系数多项式，根向量中的复数必须共轭成对。）要形成实系数多项式，根向量中的复数必须共轭成对。（2）含复数的根向量产生的多项式系数向量）含复数的根向量产生的多项式系数向量P 有可能存在截有可能存在截 断误差级的虚部。我们可用断误差级的虚部。我们可用real命令把虚数过滤掉。命令把虚数过滤掉。【例【例5.2.3】 计算有理式计算有理式 的的“商商”及及“余余”多

23、项式多项式p1=conv(1,0,2,conv(1,4,1,1); %计算分子多项式计算分子多项式p2=1 0 1 1;%注意缺项补零注意缺项补零q,r=deconv(p1,p2);cq=商多项式为商多项式为 ; cr=余多项式为余多项式为 ;disp(cq,poly2str(q,s)disp(cr,poly2str(r,s) 计算结果为：计算结果为：商多项式为商多项式为 s + 5余多项式为余多项式为 5 s2 + 4 s + 3 114232 sssss【例【例5.2.4】我们对中国】我们对中国1952年到年到1998年的国民总收入年的国民总收入GDP的的数据，用数据，用1、2、3和和4阶

24、多项式来拟合，并将结果用图形表示出来。阶多项式来拟合，并将结果用图形表示出来。GDP的数据文件为的数据文件为gdp.txtclear; AA=load(e:datagdp.txt);x=AA(:,1); y=AA(:,2);p1 = polyfit(x,y,1) % 用一阶多项式拟合用一阶多项式拟合p2 = polyfit(x,y,2) % 用二阶多项式拟合用二阶多项式拟合p3 = polyfit(x,y,4) % 用四阶多项式拟合用四阶多项式拟合p4 = polyfit(x,y,6) % 用六阶多项式拟合用六阶多项式拟合f1 = polyval(p1,x); % 计算一阶多项式的插值点计算一

25、阶多项式的插值点f2 = polyval(p2,x); % 计算二阶多项式的插值点计算二阶多项式的插值点f3 = polyval(p4,x); % 计算三阶多项式的插值点计算三阶多项式的插值点f4 = polyval(p4,x); % 计算四阶多项式的插值点计算四阶多项式的插值点subplot(2,2,1), plot(x,y,.,x,f1,-),title(一阶图一阶图)subplot(2,2,2), plot(x,y,.,x,f2,-),title(二阶图二阶图)subplot(2,2,3), plot(x,y,.,x,f3,-),title(三阶图三阶图)subplot(2,2,4),

26、plot(x,y,.,x,f4,-),title(四阶图四阶图)从结果图上看，四阶多项式和六阶多项式拟合的好。我们可以取四从结果图上看，四阶多项式和六阶多项式拟合的好。我们可以取四阶多项式，要注意的是，并不是阶数越多越好。阶多项式，要注意的是，并不是阶数越多越好。5.2.4 插值计算插值计算 插值与多项式拟合的不同之处在于，对给定的点多项式拟合不插值与多项式拟合的不同之处在于，对给定的点多项式拟合不必点点通过，而插值建立的函数是点点通过给定的数据点。必点点通过，而插值建立的函数是点点通过给定的数据点。MATLAB提供了丰富的插值方法，并配有样条插值工具箱。提供了丰富的插值方法，并配有样条插值工

27、具箱。 5.2.5 一维插值多项式一维插值多项式一维插值多项式的基本语法为：一维插值多项式的基本语法为：yi = interp1(x,y,xi,yi,method)这里这里x，y 为给定的数据，利用该数据建立一维插值多项式，为给定的数据，利用该数据建立一维插值多项式，xi为为待计算的插值点，待计算的插值点，yi为插值点的函数值。为插值点的函数值。method：为建立插值多项式的方法，可取：为建立插值多项式的方法，可取method = nearest：利用最近临方法建立插值多项式：利用最近临方法建立插值多项式method = linear： 利用线性方法建立插值多项式利用线性方法建立插值多项式m

28、ethod = spline： 利用样条函数方法建立插值多项式利用样条函数方法建立插值多项式method = pchip： 三阶厄米特方法建立逐段样条插值多项式三阶厄米特方法建立逐段样条插值多项式method = cubic： 利用三阶样条函数建立插值多项式利用三阶样条函数建立插值多项式【例【例5.2.5】对国民经济数据从】对国民经济数据从1955年以年以5年为间隔建立插值函年为间隔建立插值函数，并在求数，并在求1955年、年、1986年和预测年和预测1996年的数据。年的数据。B=ones(3,2)A=load(e:datazggdp.txt)x=A(:,1)y=A(:,2)ly=log(y

29、)x1=x(4,9,14,19,24,29,34,39,44)y1=y(4,9,14,19,24,29,34,39,44)xx=1955,1986,1996yy=interp1(x,y,xx,spline)B(:,1)=xx;B(:,2)=yy;B结果为：结果为：1955 910 1966 10201 1996 66851【例【例5.2.6】 在极坐标下插值的例子，给定椭圆上的四个点，用在极坐标下插值的例子，给定椭圆上的四个点，用样条插值并作图。样条插值并作图。theta=0:0.5:2*pi; % 产生自变量四个样点产生自变量四个样点y=-0.5 1 -0.5 -1 0.5 1 -0.5 0

30、.5 1 0.5 -1 -0.5 1 0.5; % 给出相应的因变量的点和初值给出相应的因变量的点和初值theta2=linspace(theta(1),theta(end),100); % 参量稠密化参量稠密化yy=spline(theta,y,theta2); %求稠密点上的插值求稠密点上的插值plot(yy(1,:),yy(2,:),b,y(1,:),y(2,:),or)5.2.6 高维插值高维插值我们以二维插值为例，其语法为：我们以二维插值为例，其语法为：ZI = interp2 (X, Y, Z, XI,YI, method)其中的参数和一维的类似。其中的参数和一维的类似。【例【例5

31、.2.7】以】以peaks函数为例，首先取较少的数据。然后对其求函数为例，首先取较少的数据。然后对其求插值函数，并比较各方法的差异。插值函数，并比较各方法的差异。x,y = meshgrid(-3:1:3);z = peaks(x,y);surf(x,y,z);title(较少点的表面图较少点的表面图) xi,yi = meshgrid(-3:0.25:3);zi1 = interp2(x,y,z,xi,yi,nearest);subplot(2,2,1),surf(xi,yi,zi1),title(nearest方法方法)zi2 = interp2(x,y,z,xi,yi,bilinear)

32、;subplot(2,2,2),surf(xi,yi,zi2),title(bilinear方法方法)zi3 = interp2(x,y,z,xi,yi,cubic);subplot(2,2,3),surf(xi,yi,zi3),title(bicubict方法方法)zi4 = interp2(x,y,z,xi,yi);subplot(2,2,4),surf(xi,yi,zi4),title(系统内定方法系统内定方法)5.3 微积分微积分 在大学生数学建模竞赛中微积分知识出现的比率是很高的，也在大学生数学建模竞赛中微积分知识出现的比率是很高的，也是很基础知识。比如是很基础知识。比如2002年大

33、学生数模竞赛中，年大学生数模竞赛中，A题就是一个题就是一个有关微积分知识的应用题，即有关微积分知识的应用题，即车灯线光源的优化设计车灯线光源的优化设计的问的问题。我们知道车灯是一个抛物线型结构，光源在旋转抛物面的题。我们知道车灯是一个抛物线型结构，光源在旋转抛物面的焦点处且形状为小长条形。我们要研究光源条的长短对灯光照焦点处且形状为小长条形。我们要研究光源条的长短对灯光照射的影响等等。射的影响等等。5.3.1 函数的数值导数与切平面函数的数值导数与切平面函数曲面的法线函数曲面的法线 用用i，j，k 表示直角坐标系的单位向量，曲面表示直角坐标系的单位向量，曲面在（在（x0, y0）处的法线（）处

34、的法线（Normal）n 是是在在 MATLAB 中计算与绘制法线的命令为：中计算与绘制法线的命令为：NX，NY，NZ = surfnorm（X，Y，Z） 122 yfxfkjyfixfn在在 MATLAB 中计算与绘制法线的命令为：中计算与绘制法线的命令为：NX，NY，NZ = surfnorm（X，Y，Z） 【例【例5.3.1】曲面法线演示。】曲面法线演示。y=-1:0.1:1;x=2*cos(asin(y); % 旋转曲面的旋转曲面的“母线母线”X,Y,Z=cylinder(x,20); % 形成旋转曲面形成旋转曲面surfnorm(X(:,11:21),Y(:,11:21),Z(:,1

35、1:21); %在曲面上画法线在曲面上画法线view(120,18) %控制观察角控制观察角 5.3.2 偏导数和梯度偏导数和梯度 函数函数 的偏导数定义如下：的偏导数定义如下：),(yxfxyxfyxxfyxfxx ),(),(lim),(0全微分可借助梯度（全微分可借助梯度（Gradient）定义为：）定义为： dydxfdyyfdxxfyxdf),(数值差分命令的基本语法：数值差分命令的基本语法：DX=diff（X） 求求 X 相邻元素的一阶差分相邻元素的一阶差分DX=diff（X，n） 求求X相邻元素的相邻元素的 n 阶差分阶差分DX=diff（X，n，dim）命令维）命令维 dim

36、上上, 求相邻元素的求相邻元素的n阶差分阶差分差分的含义为：差分的含义为：设序列设序列则一阶差分的定义为：则一阶差分的定义为：n阶差分依次类推。阶差分依次类推。 nxxxX,21 12312, nnxxxxxxX【例【例5.3.2】差分方法应用十分广泛，例如在经济数据中一阶差分表示发展速度，】差分方法应用十分广泛，例如在经济数据中一阶差分表示发展速度，二阶差分表示经济发展的加速度。对我国二阶差分表示经济发展的加速度。对我国5298年国民经济数据年国民经济数据GDP发展趋势发展趋势进行分析，求其一、二阶差分。并作图。进行分析，求其一、二阶差分。并作图。A=load(e:datagdp.txt)B

37、=A(:,2)B1=diff(B);B2=diff(B,2)subplot(1,3,1),plot(B),title(发展趋势发展趋势)subplot(1,3,2),plot(B1),title(发展速度发展速度)subplot(1,3,3),plot(B2),title(发展加速度发展加速度)梯度梯度:基本语法：基本语法：FX，FY=gradient（F，h） 求二元函数的梯度求二元函数的梯度FX，FY，FZ， = gradient（F，h1, h2, ） 【例【例5.3.3】 对函数对函数 作梯度图作梯度图 v = -2:0.2:2; x,y = meshgrid(v);z = x .*

38、exp(-x.2 - y.2);subplot(1,2,1);surf(x,y,z),subplot(1,2,2); px,py = gradient(z,.2,.2);contour(x,y,z), hold on, quiver(x,y,px,py) 22yxxe 方向法线的绘图命令：基本语法为：方向法线的绘图命令：基本语法为：quiver3(Z,U,V,W)quiver3(X,Y,Z,U,V,W)quiver3(.,scale)quiver3(.,LineSpec)这里：这里：X，Y，Z： 为曲线上的点。为曲线上的点。 U，V，W： 为每一点上切面的法线方向。为每一点上切面的法线方向。

39、scale： 为实数值，可正可负。表示方向箭头的指向。为实数值，可正可负。表示方向箭头的指向。【例【例5.3.4】求函数的曲面图，且作出每点的切面法向量图形。】求函数的曲面图，且作出每点的切面法向量图形。X,Y = meshgrid(-2:0.25:2,-1:0.2:1); % 产生网格点产生网格点Z = X.* exp(-X.2 - Y.2); % 计算每个网格点的函数值计算每个网格点的函数值U,V,W = surfnorm(X,Y,Z); % 计算每个点的法向量分量计算每个点的法向量分量quiver3(X,Y,Z,U,V,W,0.5); % 作法向量的图形作法向量的图形hold onsur

40、f(X,Y,Z);colormap hsvview(-35,45)axis (-2 2 -1 1 -.6 .6)hold off【例【例5.3.5】综合应用。我们来作一个旋转抛物面上某点的切面和】综合应用。我们来作一个旋转抛物面上某点的切面和该切面的法线。这里用到了作方向法线的命令。该切面的法线。这里用到了作方向法线的命令。x=0:0.1:1;y=sqrt(2*x); % 旋转曲面的旋转曲面的“母线母线”X,Y,Z=cylinder(y,20); % 形成旋转曲面形成旋转曲面U,V,W = surfnorm(X,Y,Z);Z1=Z(4)-(U(4)*(X-X(4)+V(4)*(Y-Y(4)/W

41、(4) % 求点求点X(4),Y(4),Z(4)的切平面的切平面surf(X,Y,Z); % 在曲面上画法线在曲面上画法线hold on;surf(X,Y,Z1)；hold on 在上面点处画切平面在上面点处画切平面plot3(X(4),Y(4),Z(4),r.,MarkerSize,30) % 将该点突出放大并用红色表示将该点突出放大并用红色表示U(1:3)=NaN;U(5:end)=NaNhold on;quiver3(X,Y,Z,U,V,W,-4) % 作该点处的法线，长度为作该点处的法线，长度为4,反方向反方向axis equal % 该命令是必须的，使得视觉正确该命令是必须的，使得视

42、觉正确5.3.3 数值积分数值积分 对于形如的积分对于形如的积分积分的一般命令为：积分的一般命令为：q = quad(fun,a,b)q = quad(fun,a,b,tol)q = quad(fun,a,b,tol,trace)q = quad(fun,a,b,tol,trace,p1,p2,.)q,fcnt = quadl(fun,a,b,.)其中其中fun表示被积函数，系数表示被积函数，系数a，b为积分的上下限，为积分的上下限，tol为用户给为用户给定的误差限，如不给出则系统内定的值为定的误差限，如不给出则系统内定的值为0.0000001。【例【例5.3.6】以下是一个比较著名的例子，计

43、算积分】以下是一个比较著名的例子，计算积分Q = quad(1./(1+25*x.2),-1,1) ans = 0.5494 1122511dxx也可以使用内联函数定义被积函数也可以使用内联函数定义被积函数F = inline(1./(1+25*x.2); % 定义一个内联函数定义一个内联函数Q = quad(F,-1，1); ans = 0.5494【例【例5.3.7】对空间曲线求积分，设空间曲线的参数方程为：】对空间曲线求积分，设空间曲线的参数方程为：x = t.*sin(t), y = t.*cos(t), z=t t = 0:pi/50:10*pi;plot3(t.*sin(t),t.

44、*cos(t),t,LineWidth,2)f = inline(sqrt(4*(t.*sin(t).2 + (t.*cos(t).2 + 1);len = quad(f,0,pi)计算结果为：计算结果为：len = 8.51915.3.4 二重积分二重积分二重积分的原理与一重积分一样，使用也很方便。其语法为：二重积分的原理与一重积分一样，使用也很方便。其语法为：q = dblquad(fun,xmin,xmax,ymin,ymax)q = dblquad(fun,xmin,xmax,ymin,ymax,tol)q = dblquad(fun,xmin,xmax,ymin,ymax,tol,m

45、ethod)【例【例3.3.8】求积分】求积分q=dblquad(exp(-x.2.-y2),-1,1,-1,1) q = 2.2310 1111)(22dxeyx5.4 求函数的根求函数的根 函数求根或称求函数的根，是初等数学就学习过的。然而对函数求根或称求函数的根，是初等数学就学习过的。然而对于复杂的函数、方程组求根的运算是用搜索的方法根据函数的于复杂的函数、方程组求根的运算是用搜索的方法根据函数的特性逐步接近某个根，在使用特性逐步接近某个根，在使用MATLAB命令时应该对函数有命令时应该对函数有所了解，给出某个根附近的初始点或根所在的某一个区间。所了解，给出某个根附近的初始点或根所在的某

46、一个区间。 主要求根命令为：主要求根命令为：roots：求多项式的根：求多项式的根fzero：求一般一元函数的根：求一般一元函数的根fsolve：求非线性方程组的根：求非线性方程组的根5.4.1 求多项式和一般一元函数的根求多项式和一般一元函数的根【例【例5.4.1】求多项式】求多项式 的根的根clear;p = 1,-6,-72,-27r = roots(p) p = 1 -6 -72 -27r = 12.1229 -5.7345 -0.3884 2772623 xxx求一元函数的求一元函数的0点命令点命令 fzero其基本语法为：其基本语法为：x = fzero(fun,x0)x = fz

47、ero(fun,x0,options)x = fzero(fun,x0,options,P1,P2,.)x,fval = fzero(.)x,fval,exitflag = fzero(.)x,fval,exitflag,output = fzero(.)这里这里fun： 被求解的函数，可用被求解的函数，可用inline或字符串定义或字符串定义x0： 为初初始值，可以是一个数，也可以是一个区间为初初始值，可以是一个数，也可以是一个区间fval： 目标函数值目标函数值options：输出显示项，可选择：输出显示项，可选择display： off : 没有输出没有输出iter： 将每一步的迭代结果

48、都输出来将每一步的迭代结果都输出来notify： 只有当迭代不收敛时才输出警告信息（系统内定值）只有当迭代不收敛时才输出警告信息（系统内定值）TolX： 迭代终止精度。迭代终止精度。【例【例5.4.2】在】在1，2范围范围sinx2的根并显示详细迭代步骤。的根并显示详细迭代步骤。options=optimset(Display,iter)x = fzero(sin(x*x),1 2,options)结果为：结果为： Func-count x f(x) Procedure 1 1 0.841471 initial 2 2 -0.756802 initial 3 1.52649 0.725271

49、interpolation 4 1.75821 0.0502804 interpolation 5 1.77439 -0.00687537 interpolation 6 1.77245 3.0226e-005 interpolation 7 1.77245 1.62826e-008 interpolation 8 1.77245 -1.2098e-015 interpolation 9 1.77245 1.89882e-015 interpolationZero found in the interval: 1, 2.x = 1.7725对于任意对于任意f(x) 可能有可能有0点，可能没有，

50、可能有多个点，可能没有，可能有多个0点，甚至可能点，甚至可能有无穷的有无穷的0点。因此找不到一个通解，我们可以通过人机对话，将点。因此找不到一个通解，我们可以通过人机对话，将函数显示在屏幕上，在可能的点上用鼠标采样得到初始点，然后让函数显示在屏幕上，在可能的点上用鼠标采样得到初始点，然后让计算机根据这些点自动分别地求出每一个根。计算机根据这些点自动分别地求出每一个根。【例【例5.4.3】编制程序，对任何函数求根，在运行该程序时将显示】编制程序，对任何函数求根，在运行该程序时将显示函数与零点的图形，用户在可能的零点附近用鼠标点击采点，程序函数与零点的图形，用户在可能的零点附近用鼠标点击采点，程序

51、通过循环计算每个采集点的根，并输出结果。通过循环计算每个采集点的根，并输出结果。以函数以函数 求零点为例，综合叙述相关命令的用法。求零点为例，综合叙述相关命令的用法。-10-50510-5-4-3-2-101ty(t) tbettfat )(sin)(2y=inline(sin(t)2*exp(-a*t)-b*abs(t),t,a,b); % 建立内联函数建立内联函数a=0.1;b=0.5;t=-10:0.01:10; % 对自变量采样，采样步长不宜太大。对自变量采样，采样步长不宜太大。y_char=vectorize(y); Y=feval(y_char,t,a,b);clf% 作人机对话界

52、面作人机对话界面plot(t,Y,r);hold on,plot(t,zeros(size(t),k); %画坐标横轴画坐标横轴xlabel(t);ylabel(y(t),hold off title(用鼠标在函数零点附近连续点五点用鼠标在函数零点附近连续点五点)% 循环使用求解命令循环使用求解命令fzero，求鼠标采集点附近的函数根。，求鼠标采集点附近的函数根。for i=1:5t1,y1=ginput(1);ttt(i),yyy(i),exitflag=fzero(y,t1,a,b) % 在在ttt（i）附近搜索）附近搜索0点点endA(:,1)=ttt;A(:,2)=yyy;A 则计算结

53、果为：则计算结果为：A = x f（x） -2.0074 0.0000 -2.0074 0.0000 -0.5198 0.0000 1.6738 0.0000 1.6738 -0.00005.4.2 求多元非线性方程组的根求多元非线性方程组的根求多元非线性方程组的根，算法较为复杂，我们可以利用求多元非线性方程组的根，算法较为复杂，我们可以利用fsolve命令来求解，其命令语法为：命令来求解，其命令语法为：x = fsolve(fun,x0)x = fsolve(fun,x0,options)x = fsolve(fun,x0,options,P1,P2, . )x,fval = fsolve(

54、.)x,fval,exitflag = fsolve(.)x,fval,exitflag,output = fsolve(.)x,fval,exitflag,output,jacobian = fsolve(.)这里这里fun：多元非线性方程组：多元非线性方程组x0： 一向量，初始值一向量，初始值fval： 目标函数值目标函数值output: 输出内容控制输出内容控制iterations：迭代次数输出：迭代次数输出funcCount： 计算函数值计算函数值algorithm： 算法使用算法使用cgiterations： PCG迭代次数迭代次数（仅对大型问题）（仅对大型问题）stepsize：

55、最终步长的选择最终步长的选择（仅对中型问题）（仅对中型问题）firstorderopt options： 输出显示项，可选择输出显示项，可选择Diagnostics： 输出诊断信息输出诊断信息display： off : 没有输出没有输出iter： 将每一步的迭代结果都输出来将每一步的迭代结果都输出来notify： 只有当迭代不收敛时才输出警告信息（系统内定值）只有当迭代不收敛时才输出警告信息（系统内定值）Jacobian: 雅克比方法雅克比方法on : 精确雅克比方法精确雅克比方法off： 近似雅克比方法近似雅克比方法MaxFunEvals：函数的最大上限：函数的最大上限MaxIter： 最

56、大迭代次数最大迭代次数TolFun： 终止容忍函数值终止容忍函数值TolX： 终止迭代容忍度终止迭代容忍度exitflag： 退出条件退出条件 0 表示函数收敛到某解，正常退出表示函数收敛到某解，正常退出 达到最大给定迭代次数，退出达到最大给定迭代次数，退出 0 求解过程不收敛求解过程不收敛 【例【例5.4.4】求解二元函数方程组】求解二元函数方程组 在在-5，-5为初始值的零点。为初始值的零点。F = 2*x(1) - x(2) - exp(-x(1), -x(1) + 2*x(2) - exp(-x(2);x0 = -5; -5; % 初始值初始值options=optimset(Disp

57、lay,iter); % 为为Option 设置值设置值x,fval = fsolve(F,x0,options) % 求零点求零点 02,02),(21212211xxexxyxfexxyxfx = 0.5671 0.5671fval = 1.0e-008 * -0.5319 -0.53195.5求函数极值及最优化问题求函数极值及最优化问题 MATLAB提供一些求函数极值的函数，如单变量函数求极小提供一些求函数极值的函数，如单变量函数求极小fmins和多变量函数求极小和多变量函数求极小fminsearch等。这两个函数的语法等。这两个函数的语法和参数与求根函数类似不再缀述。我们也可以打开最优

58、化工具箱使和参数与求根函数类似不再缀述。我们也可以打开最优化工具箱使用里面的更专业的方法。现将主要的命令列表如下：用里面的更专业的方法。现将主要的命令列表如下：5.5.1 无约束最优化问题无约束最优化问题【例【例5.5.1】求的】求的 最小值，给定初始值最小值，给定初始值-1x,fval,exitflag=fminbnd(x2,-1,x,optimset(TolX,1e-12,Display,off)结果为：结果为： x = 0fval = 0exitflag = 12x【例【例5.5.2】一个典型的函数】一个典型的函数是是Rosenbrock banana 函数，其典型的最小值坐标我们已经函

59、数，其典型的最小值坐标我们已经知道为（知道为（1，1）。首先将它的图形作出。）。首先将它的图形作出。（1）作）作Rosenbrock banana 函数的三维等高图函数的三维等高图x=-3:0.1:3;y=-2:0.1:4; % -3，3-2，4区域取点区域取点 X,Y=meshgrid(x,y); % 根据点在区域打上网格根据点在区域打上网格F=100*(Y-X.2).2+(1-X).2; % 每一点处计算函数值每一点处计算函数值contour3(X,Y,F,300), % 作三维等高线图形作三维等高线图形xlabel(x),ylabel(y),axis(-3,3,-2,4,0,inf),v

60、iew(161,45)hold on,plot3(1,1,0,.r,MarkerSize,20),hold off我们可以看到有一条香蕉状的谷底，事实上该函数的最小值不是对我们可以看到有一条香蕉状的谷底，事实上该函数的最小值不是对所有的方法都能找到的，因为香蕉状的谷底较平坦，因此所有的方法都能找到的，因为香蕉状的谷底较平坦，因此Rosenbrock banana 函数成了著名的搜索方法的测试函数。函数成了著名的搜索方法的测试函数。（2）建立用于搜索的）建立用于搜索的Rosenbrock banana 函数的符号形式函数的符号形式ff=inline(100*(x(2)-x(1)2)2+(1-x(