1、第第6章章 MATLAB数值计算数值计算6.1 数据处理与多项式计算数据处理与多项式计算6.2 数值微积分数值微积分6.3 离散傅立叶变换离散傅立叶变换6.4 线性方程组求解线性方程组求解6.5 非线性方程与最优化问题求解非线性方程与最优化问题求解6.6 常微分方程的数值求解常微分方程的数值求解6.7 稀疏矩阵稀疏矩阵6.1 数据处理与多项式计算数据处理与多项式计算 6.1.1 数据统计与分析数据统计与分析1. 求矩阵最大元素和最小元素求矩阵最大元素和最小元素MATLAB提供的求数据序列的最大值和最小提供的求数据序列的最大值和最小值的函数分别为值的函数分别为max和和min，两个函数的调，两个
2、函数的调用格式和操作过程类似。用格式和操作过程类似。（1）求向量的最大值和最小值）求向量的最大值和最小值 y=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，如，如果果X中包含复数元素，则按模取最大值。中包含复数元素，则按模取最大值。 y,I=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，最大，最大值的序号存入值的序号存入I，如果，如果X中包含复数元素，则按模中包含复数元素，则按模取最大值。取最大值。求向量求向量X的最小值的函数是的最小值的函数是min(X)，用法和，用法和max(X)完全相同。完全相同。例例 求向量求向量x的最大值。的最大值。命令如下：命令如下：
3、x=-43,72,9,16,23,47;y=max(x) %求向量求向量x中的最大值中的最大值y,l=max(x) %求向量求向量x中的最大值及其该元素的位置中的最大值及其该元素的位置（2）求矩阵的最大值和最小值）求矩阵的最大值和最小值求矩阵求矩阵A的最大值的函数有的最大值的函数有3种调用格式，分种调用格式，分别是：别是： max(A)：返回一个行向量，向量的第：返回一个行向量，向量的第i个元个元素是矩阵素是矩阵A的第的第i列上的最大值。列上的最大值。 Y,U=max(A)：返回行向量：返回行向量Y和和U，Y向量向量记录记录A的每列的最大值，的每列的最大值，U向量记录每列最向量记录每列最大值的
4、行号。大值的行号。 max(A,dim)：dim取取1或或2。dim取取1时，时，该函数和该函数和max(A)完全相同；完全相同；dim取取2时，该时，该函数返回一个列向量，其第函数返回一个列向量，其第i个元素是个元素是A矩矩阵的第阵的第i行上的最大值。行上的最大值。 求最小值的函数是求最小值的函数是min，其用法和，其用法和max完全完全相同。相同。例例6.1 分别矩阵分别矩阵A中各列和各行元素中的最大中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。值，并求整个矩阵的最大值和最小值。（3）两个向量或矩阵对应元素的比较）两个向量或矩阵对应元素的比较 函数函数max和和min还能对两个
5、同型的向量或矩阵进行比较，还能对两个同型的向量或矩阵进行比较，调用格式为：调用格式为： U=max(A,B)：A,B是两个同型的向量或矩阵，结果是两个同型的向量或矩阵，结果U是与是与A,B同型的向量或矩阵，同型的向量或矩阵，U的每个元素等于的每个元素等于A,B对应元素的对应元素的较大者。较大者。 U=max(A,n)：n是一个标量，结果是一个标量，结果U是与是与A同型的向量或同型的向量或矩阵，矩阵，U的每个元素等于的每个元素等于A对应元素和对应元素和n中的较大者。中的较大者。min函数的用法和函数的用法和max完全相同。完全相同。 例例 求两个求两个23矩阵矩阵x, y所有同一位置上的较大元素
6、构成的所有同一位置上的较大元素构成的新矩阵新矩阵p。2. 求矩阵的平均值和中值求矩阵的平均值和中值求数据序列平均值的函数是求数据序列平均值的函数是mean，求数据序列中值的函数，求数据序列中值的函数是是median。两个函数的调用格式为：。两个函数的调用格式为：mean(X)：返回向量：返回向量X的算术平均值。的算术平均值。median(X)：返回向量：返回向量X的中值。的中值。mean(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的算术列的算术平均值。平均值。median(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的中
7、列的中值。值。mean(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于mean(A)；当；当dim为为2时，返回一个列向量，其第时，返回一个列向量，其第i个元素是个元素是A的第的第i行的算术行的算术平均值。平均值。median(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于median(A)；当；当dim为为2时，返回一个列向量，其第时，返回一个列向量，其第i个元素是个元素是A的第的第i行的行的中值。中值。3. 矩阵元素求和与求积矩阵元素求和与求积数据序列求和与求积的函数是数据序列求和与求积的函数是sum和和prod，其使用方法类似。设其使用方法类似。设X是
8、一个向量，是一个向量，A是一是一个矩阵，函数的调用格式为：个矩阵，函数的调用格式为：sum(X)：返回向量：返回向量X各元素的和。各元素的和。prod(X)：返回向量：返回向量X各元素的乘积。各元素的乘积。sum(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的元素和。列的元素和。prod(A)：返回一个行向量，其第：返回一个行向量，其第i个元素是个元素是A的第的第i列的元素乘积。列的元素乘积。sum(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于sum(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素
9、是A的第的第i行的各元素之和。行的各元素之和。prod(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于prod(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素是A的第的第i行的各元素乘积。行的各元素乘积。例例6.2 求矩阵求矩阵A的每行元素的乘积和全部元素的每行元素的乘积和全部元素的乘积。的乘积。4. 矩阵元素累加和与累乘积矩阵元素累加和与累乘积在在MATLAB中，使用中，使用cumsum和和cumprod函数能方便地求得函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式向量和矩阵元素的累加和与累乘积向量，函数的调用格
10、式为：为：cumsum(X)：返回向量：返回向量X累加和向量。累加和向量。cumprod(X)：返回向量：返回向量X累乘积向量。累乘积向量。cumsum(A)：返回一个矩阵，其第：返回一个矩阵，其第i列是列是A的第的第i列的累加和向列的累加和向量。量。cumprod(A)：返回一个矩阵，其第：返回一个矩阵，其第i列是列是A的第的第i列的累乘积列的累乘积向量。向量。cumsum(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于cumsum(A)；当当dim为为2时，返回一个矩阵，其第时，返回一个矩阵，其第i行是行是A的第的第i行的累加行的累加和向量。和向量。cumprod(A,d
11、im)：当：当dim为为1时，该函数等同于时，该函数等同于cumprod(A)；当当dim为为2时，返回一个向量，其第时，返回一个向量，其第i行是行是A的第的第i行的累乘行的累乘积向量。积向量。5求标准方差求标准方差在在MATLAB中，提供了计算数据序列的标准方差的函数中，提供了计算数据序列的标准方差的函数std。对于向量对于向量X，std(X)返回一个标准方差。对于矩阵返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵返回一个行向量，它的各个元素便是矩阵A各列或各列或各行的标准方差。各行的标准方差。std函数的一般调用格式为：函数的一般调用格式为：Y=std(A,
12、flag,dim)其中其中dim取取1或或2。当。当dim=1时，求各列元素的标准方差；当时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。时，则求各行元素的标准方差。flag取取0或或1，当，当flag=0时，按时，按1所列公式计算标准方差，当所列公式计算标准方差，当flag=1时，按时，按2所列公式计算标准方差。缺省所列公式计算标准方差。缺省flag=0，dim=1。例例6.4 对二维矩阵对二维矩阵x，从不同维方向求出其标准方差。，从不同维方向求出其标准方差。6相关系数相关系数MATLAB提供了提供了corrcoef函数，可以求出数函数，可以求出数据的相关系数矩阵。据的相关
13、系数矩阵。corrcoef函数的调用格函数的调用格式为：式为：corrcoef(X)：返回从矩阵：返回从矩阵X形成的一个相关系形成的一个相关系数矩阵。此相关系数矩阵的大小与矩阵数矩阵。此相关系数矩阵的大小与矩阵X一一样。它把矩阵样。它把矩阵X的每列作为一个变量，然后的每列作为一个变量，然后求它们的相关系数。求它们的相关系数。corrcoef(X,Y)：在这里，：在这里，X,Y是向量，它们与是向量，它们与corrcoef(X,Y)的作用一样。的作用一样。例例6.5 生成满足正态分布的生成满足正态分布的100005随机矩随机矩阵，然后求各列元素的均值和标准方差，阵，然后求各列元素的均值和标准方差，
14、再求这再求这5列随机数据的相关系数矩阵。列随机数据的相关系数矩阵。命令如下：命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)7. 排序排序MATLAB中对向量中对向量X是排序函数是是排序函数是sort(X)，函数返，函数返回一个对回一个对X中的元素按升序排列的新向量。中的元素按升序排列的新向量。sort函数也可以对矩阵函数也可以对矩阵A的各列或各行重新排序，其的各列或各行重新排序，其调用格式为：调用格式为：Y,I=sort(A,dim)其中其中dim指明对指明对A的列还是行进行排序。若的列还是行进行排序。若dim=1，则按列排；若则按列排
15、；若dim=2，则按行排。，则按行排。Y是排序后的矩是排序后的矩阵，而阵，而I记录记录Y中的元素在中的元素在A中位置。中位置。6.1.2 数据插值数据插值1. 一维数据插值一维数据插值在在MATLAB中，实现这些插值的函数是中，实现这些插值的函数是interp1，其调用格，其调用格式为：式为：Y1=interp1(X,Y,X1,method)函数根据函数根据X,Y的值，计算函数在的值，计算函数在X1处的值。处的值。X,Y是两个等长是两个等长的已知向量，分别描述采样点和样本值，的已知向量，分别描述采样点和样本值，X1是一个向量是一个向量或标量，描述欲插值的点，或标量，描述欲插值的点，Y1是一个与
16、是一个与X1等长的插值结等长的插值结果。果。method是插值方法，允许的取值有是插值方法，允许的取值有linear、nearest、cubic、spline。注意：注意：X1的取值范围不能超出的取值范围不能超出X的给定范围，否则，的给定范围，否则，会给出会给出“NaN”错误。错误。 例例6.7 给出概率积分的数据表如表给出概率积分的数据表如表6.1所示，用不同的插值所示，用不同的插值方法计算方法计算f(0.472)。 例例6.8 某检测参数某检测参数f随时间随时间t的采样结果如表的采样结果如表5.1，用数据插，用数据插值法计算值法计算t=2，7，12，17，22，17，32，37，42，47
17、，52，57时的时的f值。值。2. 二维数据插值二维数据插值在在MATLAB中，提供了解决二维插值问题的函数中，提供了解决二维插值问题的函数interp2，其调用格式为：，其调用格式为：Z1=interp2(X,Y,Z,X1,Y1,method)其中其中X,Y是两个向量，分别描述两个参数的采样点，是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的函数值，是与参数采样点对应的函数值，X1,Y1是两个向是两个向量或标量，描述欲插值的点。量或标量，描述欲插值的点。Z1是根据相应的插是根据相应的插值方法得到的插值结果。值方法得到的插值结果。 method的取值与一维插的取值与一维插值函数相同
18、。值函数相同。X,Y,Z也可以是矩阵形式。也可以是矩阵形式。同样，同样，X1,Y1的取值范围不能超出的取值范围不能超出X,Y的给定范围，的给定范围，否则，会给出否则，会给出“NaN”错误。错误。例例6.9 设设z=x2+y2，对，对z函数在函数在0,10,2区域内进行区域内进行插值。插值。例例6.10 某实验对一根长某实验对一根长10米的钢轨进行热源的温度米的钢轨进行热源的温度传播测试。用传播测试。用x表示测量点表示测量点0:2.5:10(米米)，用，用h表示表示测量时间测量时间0:30:60(秒秒)，用，用T表示测试所得各点的温表示测试所得各点的温度度()。试用线性插值求出在一分钟内每隔。试
19、用线性插值求出在一分钟内每隔10秒、秒、钢轨每隔钢轨每隔0.5米处的温度。米处的温度。6.1.3 曲线拟合曲线拟合在在MATLAB中，用中，用polyfit函数来求得最小二乘拟合多项式的函数来求得最小二乘拟合多项式的系数，再用系数，再用polyval函数按所得的多项式计算所给出的点上函数按所得的多项式计算所给出的点上的函数近似值。的函数近似值。polyfit函数的调用格式为：函数的调用格式为：P,S=polyfit(X,Y,m)函数根据采样点函数根据采样点X和采样点函数值和采样点函数值Y，产生一个，产生一个m次多项式次多项式P及其在采样点的误差向量及其在采样点的误差向量S。其中。其中X,Y是两
20、个等长的向量，是两个等长的向量，P是一个长度为是一个长度为m+1的向量，的向量，P的元素为多项式系数。的元素为多项式系数。polyval函数的功能是按多项式的系数计算函数的功能是按多项式的系数计算x点多项式的值。点多项式的值。例例6.11 用一个用一个3次多项式在区间次多项式在区间0，2内逼近函数。内逼近函数。命令如下：命令如下：X=linspace(0,2*pi,50);Y=sin(X);P=polyfit(X,Y,3) %得到得到3次多项式的系数和误差次多项式的系数和误差6.1.4 多项式计算多项式计算1. 多项式的四则运算多项式的四则运算（1）多项式的加减运算）多项式的加减运算（2）多项
21、式乘法运算）多项式乘法运算函数函数conv(P1,P2)用于求多项式用于求多项式P1和和P2的乘积。的乘积。这里，这里，P1、P2是两个多项式系数向量。是两个多项式系数向量。（3）多项式除法）多项式除法函数函数Q,r=deconv(P1,P2)用于对多项式用于对多项式P1和和P2作除作除法运算。其中法运算。其中Q返回多项式返回多项式P1除以除以P2的商式，的商式，r返返回回P1除以除以P2的余式。这里，的余式。这里，Q和和r仍是多项式系数仍是多项式系数向量。向量。deconv是是conv的逆函数，即有的逆函数，即有P1=conv(P2,Q)+r。2. 多项式的导函数多项式的导函数对多项式求导数
22、的函数是：对多项式求导数的函数是：p=polyder(P)：求多项式：求多项式P的导函数的导函数p=polyder(P,Q)：求：求PQ的导函数的导函数p,q=polyder(P,Q)：求：求P/Q的导函数，导函数的分的导函数，导函数的分子存入子存入p，分母存入，分母存入q。上述函数中，参数上述函数中，参数P,Q是多项式的向量表示，结果是多项式的向量表示，结果p,q也是多项式的向量表示。也是多项式的向量表示。3. 多项式求值多项式求值MATLAB提供了两种求多项式值的函数：提供了两种求多项式值的函数：polyval与与polyvalm，它们的输入参数均为多项式系数向量，它们的输入参数均为多项式
23、系数向量P和自变量和自变量x。两者的区别在于前者是代数多项式求。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。值，而后者是矩阵多项式求值。(1)代数多项式求值代数多项式求值polyval函数用来求代数多项式的值，其调用函数用来求代数多项式的值，其调用格式为：格式为：Y=polyval(P,x)若若x为一数值，则求多项式在该点的值；若为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。元素求其多项式的值。例例6.14 已知多项式已知多项式x4+8x3-10，分别取，分别取x=1.2和和一个一个23矩阵为自变
24、量计算该多项式的值。矩阵为自变量计算该多项式的值。(2)矩阵多项式求值矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格式与函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。相同，但含义不同。polyvalm函数要求函数要求x为方阵，它以方为方阵，它以方阵为自变量求多项式的值。设阵为自变量求多项式的值。设A为方阵，为方阵，P代表多项式代表多项式x3-5x2+8，那么，那么polyvalm(P,A)的含义是：的含义是：A*A*A-5*A*A+8*eye(size(A)而而polyval(P,A)的含义是：的含义是：A.*A.*A-5*A.*A+8*ones(si
25、ze(A)例例6.15 仍以多项式仍以多项式x4+8x3-10为例，取一个为例，取一个22矩阵为自变量矩阵为自变量分别用分别用polyval和和polyvalm计算该多项式的值。计算该多项式的值。4. 多项式求根多项式求根n次多项式具有次多项式具有n个根，当然这些根可能是实个根，当然这些根可能是实根，也可能含有若干对共轭复根。根，也可能含有若干对共轭复根。MATLAB提供的提供的roots函数用于求多项式的函数用于求多项式的全部根，其调用格式为：全部根，其调用格式为：x=roots(P)其中其中P为多项式的系数向量，求得的根赋给向为多项式的系数向量，求得的根赋给向量量x，即，即x(1),x(2
26、),x(n)分别代表多项式的分别代表多项式的n个根。个根。例例6.16 求多项式求多项式x4+8x3-10的根。的根。命令如下：命令如下：A=1,8,0,0,-10;x=roots(A)若已知多项式的全部根，则可以用若已知多项式的全部根，则可以用poly函数建立起函数建立起该多项式，其调用格式为：该多项式，其调用格式为：P=poly(x)若若x为具有为具有n个元素的向量，则个元素的向量，则poly(x)建立以建立以x为其为其根的多项式，且将该多项式的系数赋给向量根的多项式，且将该多项式的系数赋给向量P。例例6.17 已知已知 f(x)(1) 计算计算f(x)=0 的全部根。的全部根。(2) 由
27、方程由方程f(x)=0的根构造一个多项式的根构造一个多项式g(x)，并，并与与f(x)进行对比。进行对比。命令如下：命令如下：P=3,0,4,-5,-7.2,5;X=roots(P) %求方程求方程f(x)=0的根的根G=poly(X) %求多项式求多项式g(x)6.2 数值微积分数值微积分6.2.1 数值微分数值微分1. 数值差分与差商数值差分与差商2. 数值微分的实现数值微分的实现在在MATLAB中，没有直接提供求数值导数的函数，只有计中，没有直接提供求数值导数的函数，只有计算向前差分的函数算向前差分的函数diff，其调用格式为：，其调用格式为：DX=diff(X)：计算向量：计算向量X的
28、向前差分，的向前差分，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，按行计算差分。，按行计算差分。例例6.18 设设x由由0,2间均匀分布的间均匀分布的10个点组成，求个点组成，求sinx的的13阶差分。阶差分。命令如下：命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(
29、Y); %计算计算Y的一阶差分的一阶差分D2Y=diff(Y,2); %计算计算Y的二阶差分的二阶差分，也可用命令，也可用命令diff(DY)计算计算D3Y=diff(Y,3); %计算计算Y的三阶差分的三阶差分，也可用，也可用diff(D2Y)或或diff(DY,2)例例6.19 用不同的方法求函数用不同的方法求函数f(x)的数值导数，并在同一个坐的数值导数，并在同一个坐标系中做出标系中做出f(x)的图像。的图像。程序如下：程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3
30、+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %用用5次多项式次多项式p拟合拟合f(x)dp=polyder(p); %对拟合多项式对拟合多项式p求导数求导数dpdpx=polyval(dp,x); %求求dp在假设点的函数值在假设点的函数值dx=diff(f(x,3.01)/0.01; %直接对直接对f(x)求数值导数求数值导数gx=g(x); %求函数求函数f的导函数的导函数g在假设点的导数在假设点的导数plot(x,dpx,x,dx,.,x,gx,-); %作图作图6.2.2 数值积分数值积分1.
31、数值积分基本原理数值积分基本原理 求解定积分的数值方法多种多样，如简单求解定积分的数值方法多种多样，如简单的梯形法、辛普生的梯形法、辛普生(Simpson) 法、牛顿法、牛顿柯特斯柯特斯(Newton-Cotes)法等都是经常采用的法等都是经常采用的方法。它们的基本思想都是将整个积分区方法。它们的基本思想都是将整个积分区间间a,b分成分成n个子区间个子区间xi,xi+1，i=1,2,n，其中其中x1=a，xn+1=b。这样求定积分问题就分。这样求定积分问题就分解为求和问题。解为求和问题。2. 数值积分的实现数值积分的实现(1)被积函数是一个解析式被积函数是一个解析式MATLAB提供了提供了qu
32、ad函数和函数和quadl函数来求函数来求定积分。它们的调用格式为：定积分。它们的调用格式为： quad(filename,a,b,tol,trace) quadl(filename,a,b,tol,trace)例例6.20 用两种不同的方法求定积分。用两种不同的方法求定积分。先建立一个函数文件先建立一个函数文件ex.m：function ex=ex(x)ex=exp(-x.2);然后在然后在MATLAB命令窗口命令窗口，输入命令：，输入命令：format longI=quad(ex,0,1) %注意函数名应加字符引号注意函数名应加字符引号I = 0.74682418072642I=quadl
33、(ex,0,1)I = 0.74682413398845也可不建立关于被积函数的函数文件，而使用语句函数也可不建立关于被积函数的函数文件，而使用语句函数(内联函数内联函数)求解，命令如求解，命令如下：下：g=inline(exp(-x.2); %定义一个语句函数定义一个语句函数g(x)=exp(-x2)I=quadl(g,0,1) %注意函数名不加注意函数名不加号号I = 0.74682413398845format short(2)被积函数由一个表格定义被积函数由一个表格定义在科学实验和工程应用中，函数关系往往是不知道在科学实验和工程应用中，函数关系往往是不知道的，只有实验测定的一组样本点和
34、样本值，这时，的，只有实验测定的一组样本点和样本值，这时，就无法使用就无法使用quad函数计算其定积分。在函数计算其定积分。在MATLAB中，对由表格形式定义的函数关系的求定积分问中，对由表格形式定义的函数关系的求定积分问题用题用trapz(X,Y)函数。其中向量函数。其中向量X、Y定义函数关定义函数关系系Y=f(X)。X、Y是两个等长的向量：是两个等长的向量：X=(x1,x2,xn)，Y=(y1,y2,yn)，并且，并且x1x2 cholMatrix must be positive definite命令执行时，出现错误信息，说明命令执行时，出现错误信息，说明A为非正定矩阵。为非正定矩阵。6
35、.4.2 迭代解法迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括中，迭代解法主要包括 Jacobi迭代法、迭代法、Gauss-Serdel迭代迭代法、超松弛迭代法和两步迭代法。法、超松弛迭代法和两步迭代法。1Jacobi迭代法迭代法对于线性方程组对于线性方程组Ax=b，如果，如果A为非奇异方阵，即为非奇异方阵，即aii0(i=1,2,n)，则可将，则可将A分解为分解为A=D-L-U，其中，其中D为对为对角阵，其元素为角阵，其元素为A的对角元素，的对角元素，L与与U为为A的下三角阵和上的下三角阵和上三角阵，于是
36、三角阵，于是Ax=b化为：化为：x=D-1(L+U)x+D-1b与之对应的迭代公式为：与之对应的迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b这就是这就是Jacobi迭代公式。如果序列迭代公式。如果序列x(k+1)收敛于收敛于x，则，则x必必是方程是方程Ax=b的解。的解。Jacobi迭代法的迭代法的MATLAB函数文件函数文件Jacobi.m如下：如下：function y,n=jacobi(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y; y=B*x0+f; n=n+1;end例例6.28 用用Jacobi迭代
37、法求解线性方程组。设迭代初值为迭代法求解线性方程组。设迭代初值为0，迭代精度为迭代精度为10-6。在命令中调用函数文件在命令中调用函数文件Jacobi.m，命令如下：，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=jacobi(A,b,0,0,0,1.0e-6)2Gauss-Serdel迭代法迭代法在在Jacobi迭代过程中，计算时，已经得到，不必再用，即原迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：，于是得到：x
38、(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为该式即为Gauss-Serdel迭代公式。和迭代公式。和Jacobi迭代相比，迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。迭代用新分量代替旧分量，精度会高些。Gauss-Serdel迭代法的迭代法的MATLAB函数文件函数文件gauseidel.m如下：如下：function y,n=gauseidel(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y; y=G*x0+f; n=n+1;end例例6.29 用用Gauss-Serdel迭代法求解
39、下列线性方程组。设迭代迭代法求解下列线性方程组。设迭代初值为初值为0，迭代精度为，迭代精度为10-6。在命令中调用函数文件在命令中调用函数文件gauseidel.m，命令如下：，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=gauseidel(A,b,0,0,0,1.0e-6)例例6.30 分别用分别用Jacobi迭代和迭代和Gauss-Serdel迭代法求解下列线迭代法求解下列线性方程组，看是否收敛。性方程组，看是否收敛。命令如下：命令如下：a=1,2,-2;1,1,1;2,2,1;b=9;7;6;x,n=jacobi(a,b,0;0;0)x,n=g
40、auseidel(a,b,0;0;0)6.4.3 求线性方程组的通解求线性方程组的通解线性方程组的求解分为两类：一类是求方程组的惟一解即特线性方程组的求解分为两类：一类是求方程组的惟一解即特解，另一类是求方程组的无穷解即通解。这里对线性方程解，另一类是求方程组的无穷解即通解。这里对线性方程组组 Ax=b的求解理论作一个归纳。的求解理论作一个归纳。(1)当系数矩阵当系数矩阵A是一个满秩方阵时，方程是一个满秩方阵时，方程Ax=b称为恰定方程，称为恰定方程，方程有惟一解方程有惟一解x=A-1b，这是最基本的一种情况。一般用，这是最基本的一种情况。一般用x=Ab求解速度更快。求解速度更快。(2)当方程
41、组右端向量当方程组右端向量b=0时，方程称为齐次方程组。齐次方时，方程称为齐次方程组。齐次方程组总有零解，因此称解程组总有零解，因此称解x=0为平凡解。当系数矩阵为平凡解。当系数矩阵A的的秩小于秩小于n(n为方程组中未知变量的个数为方程组中未知变量的个数)时，齐次方程组有时，齐次方程组有无穷多个非平凡解，其通解中包含无穷多个非平凡解，其通解中包含n-rank(A)个线性无关个线性无关的解向量，用的解向量，用MATLAB的函数的函数null(A,r)可求得基础解系。可求得基础解系。(3)当方程组右端向量当方程组右端向量b0时，系数矩阵的秩时，系数矩阵的秩rank(A)与其增广矩阵的秩与其增广矩阵
42、的秩rank(A,b)是判断其是否有解是判断其是否有解的基本条件：的基本条件：当当rank(A)=rank(A,b)=n时，方程组有惟一解：时，方程组有惟一解：x=Ab 或或 x=pinv(A)*b。当当rank(A)=rank(A,b)n时，方程组有无穷多个时，方程组有无穷多个解，其通解解，其通解=方程组的一个特解方程组的一个特解+对应的齐次方程对应的齐次方程组组Ax=0的通解。可以用的通解。可以用Ab求得方程组的一个特求得方程组的一个特解，用解，用null(A,r)求得该方程组所对应的齐次方程求得该方程组所对应的齐次方程组的基础解系，基础解系中包含组的基础解系，基础解系中包含n-rank(
43、A)个线性个线性无关的解向量。无关的解向量。当当rank(A)rank(A,b)时，方程组无解。时，方程组无解。 有了上面这些讨论，可以设计一个求解线有了上面这些讨论，可以设计一个求解线性方程组的函数文件性方程组的函数文件line_solution.m。在例。在例中可以调用中可以调用line_solution.m文件来解线性方文件来解线性方程组。程组。 6.5 非线性方程与最优化问题求解非线性方程与最优化问题求解6.5.1 非线性方程数值求解非线性方程数值求解1. 单变量非线性方程求解单变量非线性方程求解 在在MATLAB中提供了一个中提供了一个fzero函数，可以用来求单变量函数，可以用来求
44、单变量非线性方程的根。该函数的调用格式为：非线性方程的根。该函数的调用格式为： z=fzero(fname,x0,tol,trace)其中其中fname是待求根的函数文件名，是待求根的函数文件名，x0为搜索的起点。一个为搜索的起点。一个函数可能有多个根，但函数可能有多个根，但fzero函数只给出离函数只给出离x0最近的那个根。最近的那个根。tol控制结果的相对精度，缺省时取控制结果的相对精度，缺省时取tol=eps，trace 指定指定迭代信息是否在运算中显示，为迭代信息是否在运算中显示，为1时显示，为时显示，为0时不显示，时不显示，缺省时取缺省时取trace=0。例例6.33 求求f(x)在
45、在x0=-5和和x0=1作为迭代初值时的零点。作为迭代初值时的零点。先建立函数文件先建立函数文件fz.m：function f=fz(x)f=x-1/x+5;然后调用然后调用fzero函数求根。：函数求根。：fzero(fz,-5) %以以-5作为迭代初值作为迭代初值ans = -5.1926fzero(fz,1) %以以1作为迭代初值作为迭代初值ans = 0.1926 2. 非线性方程组的求解非线性方程组的求解 对于非线性方程组对于非线性方程组F(X)=0，用，用fsolve函数求其数值解。函数求其数值解。fsolve函数的调用格式为：函数的调用格式为： X=fsolve(fun,X0,o
46、ption)其中其中X为返回的解，为返回的解，fun是用于定义需求解的非线性方程组的是用于定义需求解的非线性方程组的函数文件名，函数文件名，X0是求根过程的初值，是求根过程的初值，option为最优化工具为最优化工具箱的选项设定。最优化工具箱提供了箱的选项设定。最优化工具箱提供了20多个选项，用户可多个选项，用户可以使用以使用optimset命令将它们显示出来。如果想改变其中某命令将它们显示出来。如果想改变其中某个选项，则可以调用个选项，则可以调用optimset()函数来完成。例如，函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中选项决定函数调用时中间结果的显示方式
47、，其中off为不显示，为不显示，iter表示每步都显示，表示每步都显示，final只显示只显示最终结果。最终结果。optimset(Display,off)将设定将设定Display选项为选项为off。例例6.34 求下列方程组在求下列方程组在(1，1，1)附近的解并对结果进行验附近的解并对结果进行验证。证。首先建立函数文件首先建立函数文件myfun.m。function F=myfun (X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+y+z2*exp(x);F(2)=x+y+z;F(3)=x*y*z;在给定的初值在给定的初值x0=1，y0=1，z0=1下，调用下，调用
48、fsolve函数求方程的根。函数求方程的根。X=fsolve(myfun,1,1,1,optimset(Display, off)X = 0.0224 -0.0224 -0.0000 6.5.2 无约束最优化问题求解无约束最优化问题求解在实际应用中，许多科学研究和工程计算问题都可以归结为在实际应用中，许多科学研究和工程计算问题都可以归结为一个最小化问题，如能量最小、时间最短等。一个最小化问题，如能量最小、时间最短等。MATLAB提供了提供了3个求最小值的函数，它们的调用格式为：个求最小值的函数，它们的调用格式为：(1)x,fval=fminbnd(filename,x1,x2,option)：
49、求一元函数在：求一元函数在(xl,x2)区间中的极小值点区间中的极小值点x和最小值和最小值fval。(2)x,fval=fminsearch(filename,x0,option)：基于单纯形算法：基于单纯形算法求多元函数的极小值点求多元函数的极小值点x和最小值和最小值fval。(3)x,fval=fminunc(filename,x0,option)：基于拟牛顿法求多：基于拟牛顿法求多元函数的极小值点元函数的极小值点x和最小值和最小值fval。MATLAB没有专门提供求函数最大值的函数，但只要注意没有专门提供求函数最大值的函数，但只要注意到到-f(x)在区间在区间(a,b)上的最小值就是上的
50、最小值就是f(x)在在(a,b)的最大值，的最大值，所以所以fminbnd(-f,x1,x2)返回函数返回函数f(x)在区间在区间(x1,x2)上的最大上的最大值。值。例例6.36 求函数在区间求函数在区间(-10，1)和和(1，10)上的最小值点。上的最小值点。首先建立函数文件首先建立函数文件fx.m：function f=f(x)f=x-1/x+5;上述函数文件也可用一个语句函数代替：上述函数文件也可用一个语句函数代替：f=inline(x-1/x+5)再在再在MATLAB命令窗口，输入命令：命令窗口，输入命令：fminbnd(fx,-10,-1) %求函数在求函数在(-10，-1)内的最
