1、第6章 MATLAB数值计算,6.1 数据处理与多项式计算 6.2 傅立叶分析 6.3 数值微积分 6.4 线性方程组求解 6.5 常微分方程的数值求解 6.6 非线性方程的数值求解 6.7 稀疏矩阵,6.1 数据处理与多项式计算,6.1.1 数据统计与分析 1. 求矩阵最大和最小元素 (1)求向量的最大最小元素 y=max(X) 返回向量X的最大元素存入y。 y,I=max(X) 返回向量X的最大元素存入y，最大元素的序号存入I。 (2)求矩阵的最大和最小元素 max(A) 返回一个行向量，向量的第i个元素是A矩阵的第i列上的最大元素。 Y,U=max(A) 返回两个行向量，Y向量记录A的每
2、列的最大元素，U向量记录每列最大元素的行号。 max(A,dim) dim取1或2。dim取1时，该函数和max(A)完全相同。dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大元素。,(3)两个向量或矩阵对应元素的比较 U=max(A,B) A,B是两个同型的向量或矩阵。结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。 U=max(A,n) n是一个标量。结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。 min函数的用法和max完全相同。,例6.1 求矩阵A的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 命
3、令如下： A=13,-56,78;25,63,-235;78,25,563;1,0,-1; max(A,2) %求每行最大元素 min(A,2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素%求整个矩阵的最大元素%求整个矩阵的最小元素,max(max(A) min(min(A),2. 求矩阵的平均值和中值求矩阵和向量元素的平均值的函数是mean ;求中值的函数是median。它们的调用方法和max函数完全相同。 3. 矩阵元素求和与求积矩阵和向量求和与求积的基本函数是sum和prod，其使用方法和max类似。,例6.2 求矩阵A的每行元素的乘积和全部元素的乘
4、积。 命令如下： A=1,2,3,4;5,6,7,8;9,10,11,12; S=prod(A,2) prod(S) %求A的全部元素的乘积4. 矩阵元素累加和与累乘积 MATLAB中，使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的用法和sum及prod相同 例6.3求向量X=(1！，2！，3！，10！)。 命令如下： X=cumprod(1:10),5. 标准方差 MATLAB中， 提供了计算数据序列的标准方差的函数std。 对于向量X，std(X)返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差
5、。std函数的一般调用格式为： std(A,FLAG,dim)其中dim取1或2。当dim=1时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。FLAG取0或1。,6. 元素排序 MATLAB中对向量X是排序函数是sort(X)，函数返回一个对X中的元素按升序排列的新向量。 sort函数也可以对矩阵A的各列(或行)重新排序，其调用格式为： Y,I=sort(A,dim)其中dim指明对A的列还是行进行排序，若dim=1，则按列排，若dim=2，则按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。,例6.4 对矩阵做各种排序。 命令如下： A=1,-8,5;4,12,6;1
6、3,7,-13; sort(A) %对A的每列按升序排序 -sort(-A,2) %对A的每行按降序排序 X,I=sort(A) %对A按列排序，并将每个元素所在行号送矩阵I,6.1.2 数值插值 1. 一维数值插值 interp1函数调用格式为： Y1=interp1(X,Y,X1,method) 函数根据X、Y的值，计算函数在X1处的值。X、Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。 method是插值方法，允许的取值有linear(线性插值)、nearest(最近插值)、spline(三次样条插值)、cubic（
7、三次多项式插值），缺省值是linear。,例6.5 用不同的插值方法计算sin(x)在/2点的值。 这是一个一维插值问题。 在MATLAB命令窗口，输入命令： X=0:0.2:pi;Y=sin(X); %给出X、Y interp1(X,Y,pi/2) %用缺省方法(即线性插值方法)计算sin(/2) interp1(X,Y,pi/2,nearest) %用最近方法计算sin(/2) interp1(X,Y,pi/2,linear) %用线性方法计算sin(/2) interp1(X,Y,pi/2,spline) %用三次样条方法计算sin(/2) interp1(X,Y,pi/2,cubic)
8、 %用三次多项式方法计算sin(/2)MATLAB中有一个专门的三次样条插值函数Y1=spline(X,Y,X1)，其功能及使用方法与函数Y1=interp1(X,Y,X1,spline)完全相同。,例6.6 已知检测参数f随时间t的采样结果，用数值插值法计算t=2，7，12，17，22，17，32，37，42，47，52，57时f的值。 这是一个一维数值插值问题，命令如下： T=0:5:65; X=2:5:57; F=3.2015,2.2560,879.5,1835.9,2968.8,4136.2,5237.9,6152.7,. 6725.3,6848.3,6403.5,6824.7,732
9、8.5,7857.6; F1=interp1(T,F,X) %用线性方法插值 F1=interp1(T,F,X,nearest) %用最近方法插值 F1=interp1(T,F,X,spline) %用三次样条方法插值 F1=interp1(T,F,X,cubic) %用三次多项式方法插值,2. 二维数值插值 MATLAB中，提供了解决二维插值问题的函数。其调用格式为： Z1=interp2(X,Y,Z,X1,Y1,method)其中，X、Y是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的采样变量的样本值，X1、Y1是两个向量或标量，描述欲插值的点。method的取值与一维插值函数
10、相同。,例6.7a 设Z=x2+y2，对Z函数在(0，1)(0，2)区域内进行插值。 命令如下： x=0:0.1:10;y=0:0.2:20; X,Y=meshgrid(x,y); Z=X.2+Y.2; interp2(x,y,Z,0.5,0.5) %对函数在(0.5,0.5)点进行插值 interp2(x,y,Z,0.5 0.6,0.4) %对函数在(0.5,0.4)点和(0.6，0.4)点进行插值 interp2(x,y,Z,0.5 0.6,0.4 0.5) %对函数在(0.5,0.4)点和(0.6，0.5)点进行插值 interp2(x,y,Z,0.5 0.6 ,0.4 0.5) %对函
11、数在(0.5,0.4),(0.6,0.4),(0.5,0.5)和(0.6,0.5)点进行插值,例6.7b 某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度()。试用线性插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。 程序如下： x=0:2.5:10; h=0:30:60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; xi=0:0.5:10; hi=0:10:60; temps=interp2(x,h,T,xi,hi,cubic); mes
12、h(xi,hi,temps);,3. 三维数值插值 对三维函数插值的函数是interp3，其使用方法和interp2相同。其调用格式为： W1=interp3(X,Y,Z,W,X1,Y1,Z1,method) 函数返回三维插值结果。其中，X、Y、Z是三个向量，分别描述三个参数的采样点，W是与参数采样点对应的采样变量的样本值，X1、Y1、Z1是三个向量或标量，描述欲插值的点。method是插值方法，可选，其缺省值是 line。method的取值与一、二维插值函数相同。,6.1.3 曲线拟合 MATLAB中，提供了解决使用最小二乘法进行曲线拟合的函数。调用格式为： P,S=polyfit(X,Y,
13、m)函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X、Y是两个等长的向量，P是一个长度为m+1的向量。,例6.8 用一个3次多项式在区间0，2内逼近函数sin(x)。 命令如下： X=linspace(0,2*pi,50); Y=sin(X); P,S=polyfit(X,Y,3) plot(X,Y, -*,X,polyval(P,X), :o),6.1.4 多项式计算 1. 多项式的建立 已知一个多项式的全部根X求多项式系数的函数是poly(X)，该函数返回以X为全部根的一个多项式P，当X是一个长度为m的向量时，P是一个长度为m+1的向量。 2. 多项式
14、求根 求多项式p(x)的根的函数是roots(P)，这里，P是p(x)的系数向量，该函数返回方程p(x)=0的全部根(含重根，复根)。 3. 多项式求值 求多项式p(x)在某点或某些点的函数值的函数是polyval(P,x)。若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。,例6.9 已知一个多项式，计算： (1)计算f(x)=0 的全部根。 (2)由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。 (3)计算f(5)、f(7.8)、f(9.6)、f(12.3)的值。 命令如下： P=3,0,4,-5,-7.2,5; X=roo
15、ts(P) %求方程f(x)=0的根 G=poly(X) %求多项式g(x) X0=5,7.8,9.6,12.3; f=polyval(P,X0) %求多项式f(x)在给定点的值多项式求值还有一个函数是polyvalm，其调用格式与 polyval相同，但含义不同。polyvalm函数要求x为方阵， 它以方阵为自变量求多项式的值。,4. 多项式的四则运算 (1)多项式的加减法 当两个多项式的次数不同时，要在一个较低次幂的多 项式系数向量前补0，使两个系数向量等长。 (2)多项式的乘法 函数conv(P1,P2)用于求多项式P1和P2的乘积。 (3)多项式的除法 函数Q,r=deconv(P1,
16、P2)用于对多项式P1和P2作除法 运算。其中Q返回多项式P1除以P2的商式，r返回P1除 以P2的余式。这里，Q和r仍是多项式系数向量。 deconv是conv的逆函数，即有P1=conv(P2,Q)+r。,例6.10 设有两个多项式，计算： (1)求f(x)+g(x)、f(x)-g(x)。 (2)求f(x)g(x)、f(x)/g(x)。 在MATLAB命令窗口，输入命令： f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,g; f+g1 %求f(x)+g(x) f-g1 %求f(x)-g(x) conv(f,g) %求f(x)*g(x) Q,r=deconv(f,g) %
17、求f(x)/g(x)，商式送Q，余式送r。,5. 多项式的导函数 对多项式求导数的函数是： p=polyder(P) 求多项式P的导函数 p=polyder(P,Q) 求P*Q的导函数 p,q=polyder(P,Q) 求P/Q的导函数，导函数的分子存入p，分母存入q。,例6.11 求有理分式的导数。 命令如下： P=3,5,0,-8,1,-5; Q=10,5,0,0,6,0,0,7,-1,0,-100; p,q=polyder(P,Q),6.1.5 函数的最大值与最小值 MATLAB中用于求最小值的函数是： fmin(f,a,b) 求单变量函数f(x)在区间(a,b)上的最小值点。 fmin
18、s(F,X0) 求多变量函数F(x)在估计值X0附近的最小值点。 MATLAB没有专门提供求函数最大值点的函数，但只 要注意到-f(x)在区间(a,b)上的最小值点就是f(x)在(a,b) 的最大值点，所以fmin(-f,a,b)返回函数f(x)在区间(a,b) 上的最大值。,例6.12 求函数f(x)在区间(-10，1)和(1，10)上的最小值点。 首先建立函数文件fx.m： function f=f(x) f=x-1/x+5; return 上述函数文件也可用一个语句函数代替： Ff=inline( x-1/x+5 ) 再在MATLAB命令窗口，输入命令： fmin(fx,-10,-1)
19、%求函数在区间(-10，-1)内的最小值点 fmin(ff,1,10) %求函数在区间(1，10)内的最小值点。注意函数名f不用加,例6.13 设有函数f(x,y,z)，求函数f在(0.5,0.5,0.5)附近的最小值。 建立函数文件fxyz.m： function f=f(u) x=u(1);y=u(2);z=u(3); f=x+y.2./x/4+z.2./y+2./z; return 在MALAB命令窗口，输入命令： U=fmins(fxyz,0.5,0.5,0.5) %求函数的最小值点 fxyz(U) %求函数的最小值,6.2 傅立叶分析,MATLAB中，提供了对向量(或直接对矩阵的行或
20、列)进行离散傅立叶变换的函数， 其调用格式是： Y=fft(X,n,dim) (1)当X是一个向量时，返回对X的离散傅立叶变换。 (2)当X是一个矩阵时，返回一个矩阵并送Y，其列(行)是对X的列(行)的离散傅立叶变换。,例6.14 求X=(1，0，-3，5，2)的离散傅立叶变换。 在MATLAB命令窗口，输入命令： X=1,0,-3,5,2; Y=fft(X) %对X进行变换3. 离散傅立叶变换的逆变换 MATLAB中，对向量(或直接对矩阵的行或列)进行 离散傅立叶逆变换的函数的调用方法是： Y=ifft(X,n,dim) 函数对X进行离散傅立叶逆变换。其中X、n、dim的 意义及用法和离散傅
21、立叶变换函数fft完全相同。,例6.15 对矩阵A的列向量、行向量分别进行离散傅立叶变换、并对变换结果进行逆变换。 命令如下： A=3,2,1,1;-5,1,0,1;3,2,1,5; fftA=fft(A) %求A的列向量的傅立叶变换 fftA2=fft(A,4,2) %求A的行向量的傅立叶变换 ifft(fftA) %对矩阵fftA的列向量进行傅立叶逆变换，结果应等于A ifft(fftA2,4,2) %对矩阵fftA2的行向量进行傅立叶逆变换，其结果应等于A,fftA =1.0000 5.0000 2.0000 7.0000 4.0000 + 6.9282i 0.5000 + 0.8660
22、i 0.5000 + 0.8660i -2.0000 + 3.4641i4.0000 - 6.9282i 0.5000 - 0.8660i 0.5000 - 0.8660i -2.0000 - 3.4641ifftA2 =7.0000 2.0000 - 1.0000i 1.0000 2.0000 + 1.0000i-3.0000 -5.0000 -7.0000 -5.0000 11.0000 2.0000 + 3.0000i -3.0000 2.0000 - 3.0000ians =3.0000 2.0000 1.0000 1.0000-5.0000 1.0000 0.0000 1.00003
23、.0000 2.0000 1.0000 5.0000 ans =3 2 1 1-5 1 0 13 2 1 5,6.3 数值微积分,6.3.1 数值微分 MATLAB中，没有直接提供求数值导数的函数，只有 计算向前差分的函数。 DX=diff(X) 计算向量X的向前差分，DX(i)=X(i+1)-X(i)，0in。 DX=diff(X,n) 计算X的n阶向前差分， diff(X,2)=diff(diff(X)。 DX=diff(A,n,dim) 计算矩阵A的n阶差分，dim=1时(缺 省状态)，按列计算差分，dim=2，按行计算差分。,例6.16 求向量sin(X)的13阶差分。设X由0，2间均
24、匀分布的10个点组成。 命令如下： X=linspace(0,2*pi,10); Y=sin(X); DY=diff(Y); %计算Y的一阶差分 D2Y=diff(Y,2); %计算Y的二阶差分，也可用命令diff(DY)计算 D3Y=diff(Y,3); %计算Y的三阶差分，也可用diff(D2Y)或diff(DY,2),例6.17 用不同的方法求函数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+2*x.
25、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求导数dp dpx=polyval(dp,x); %求dp在假设点的函数值 dx=diff(f(x,3.01)/0.01; %直接对f(x)求数值导数 gx=g(x); %求函数f的导函数g在假设点的导数 plot(x,dpx,x,dx,g.,x,gx,r-); %作图,6.3.2数值积分 (1)被积函数是一个解析式函数quad(f,a,b,tol,trace)用于求被积函数f(x)在a,
26、b上的定积分，tol是计算精度，缺省值是0.000001。trace非0时，画出积分图形。 注意，调用quad函数时，先要建立一个描述被积函数f(x)的函数文件或语句函数。 当被积函数f含有一个以上的变量时，quad函数的调用格式为：quad(f,a,b,tol,trace,g1,g2)其中 f,a,b,tol,trace等参数的含义同前。 数值积分函数还有一种形式quadl，其用法与quad完全相同。,数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分
27、成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,例6.18 用两种不同的方法求积分。 先建立一个函数文件ex.m： function ex=ex(x) ex=exp(-x.2); %注意应用点运算 return 然后，在MATLAB命令窗口，输入命令： quad(ex,0,1,1e-6) %注意函数名应加字符引号 quadl(ex,0,1,1e-6) %用另一函数求积分也可不建立关于被积函数的函数文件，而使用语句函数(内联函数)求解，命令如下： g=inline(exp(-x.2); %定义一个语句函数g(x)=exp(-x2) I=
28、quadl(g,0,1) %注意函数名不加号,(2)被积函数由一个表格定义MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X、Y定义函数关系Y=f(X)。例6.19 用trapz函数计算积分。 在MATLAB命令窗口，输入命令： X=0:0.01:1;Y=exp(-X.2); trapz(X,Y),(3)二重积分 例6.20 计算二重积分。 建立一个函数文件fixy.m： function f=f(x,y) f=exp(-x.2-y.2); return 建立一个命令文件ftxy1.m： for i=1:20int2(i)=quad(fixy,0,1
29、,x(i); %在二维函数fixy中以x=x(i)代入并对y积分。 end 在MATLAB命令窗口，输入命令： x=linspace(0,1,20); ftxy1 trapz(x,int2),实际上，MATLAB提供了计算二重积分的函数： dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在a,bc,d区域上的二重积分。 参数tol，trace的用法与函数quad完全相同。如果直接使用这里介绍的二重积分函数dblquad来求解 本例就非常简单，命令如下： g=inline(exp(-x.2-y.2); dblquad(g,0,1,0,1) %直接调用二重积分函数求解
30、,6.4 线性方程组求解,6.4.1 直接解法 1利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运算符“”求解：x=Ab 例6.21 用直接解法求解下列线性方程组。 命令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; x=Ab,6.4.2 利用矩阵的分解求解线性方程组 矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇异分解等。1. 矩阵的LU分解 MATLAB中，完成LU分解的函数是：
31、(1)L,U=lu(A) 将方阵A分解为交换下三角矩阵L和上三角矩阵U，使 A=LU。 (2)L,U,P=lu(A) 将方阵A分解为下三角矩阵L和上三角矩阵U，使 PA=LU。,实现LU分解后，线性方程组Ax=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。例6.22 用LU分解求解例6.21中的线性方程组。 命令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; L,U=lu(A); x=U(Lb) 或采用LU分解的第2种格式，命令如下： L,U ,P=lu(A); x=U(LP*b),2. 矩阵的QR分解对矩阵
32、X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为： Q,R=qr(X)：产生一个一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。 Q,R,E=qr(X)：产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。实现QR分解后，线性方程组Ax=b的解x=R(Qb)或x=E(R(Qb),例6.23 用QR分解求解例6.21中的线性方程组。 命令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; Q,R
33、=qr(A); x=R(Qb) 或采用QR分解的第2种格式，命令如下： Q,R,E=qr(A); x=E*(R(Qb),6.4.3线性方程组求解 1. 线性方程组解的一般讨论 解线性方程组的一般函数文件如下： function x,y=line_solution(A,b)m,n=size(A);y=;if norm(b)0 %非齐次方程组if rank(A)=rank(a,b) %方程组相容 if rank(A)=m %有唯一解x=Ab;else %方程组有无穷多个解,基础解系disp(原方程组有有无穷个解，其齐次方程组的基础解系为y，特解为x);y=null(A,r);x=Ab; endel
34、se %方程组不相容，给出最小二乘法解disp(方程组的最小二乘法解是：);x=Ab; endelse %齐次方程组if rank(A)=n %列满秩x=zero(m,1) %0解else %非0解disp(方程组有无穷个解，基础解系为x);x=null(A,r);endend return,2. 应用举例 例6.24 求线性方程组的解。 在MATLAB命令窗口，输入命令： A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6; x,y=line_solution(A,b) %调用自定义函数 例6.25 求下列线性方程组的解。 在MATLAB命令窗口
35、，输入命令： A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2; x,y=line_solution(A,b),6.5 常微分方程的数值求解,基于龙格库塔法，MATLAB提供了求常微分方程数 值解的函数，一般调用格式为：X,Y=ode23(f,x0,xn,y0)X,Y=ode45(f,x0,xn,y0)其中，X、Y是两个向量，X对应自变量x在求解区间x1，xn的一组采样点，其采样密度是自适应的，无需指定；Y是与X对应的一组解， f是一个函数，x0，xn代表自变量的求解区间，y0=y(x0)，由方程的初值给定。函数在求解区间x0，xn内，自动设立采样点向量X，并求出解函数y在采
36、样点X处的样本值。,例6.26 求微分方程初值问题在1，3区间内的数值解，并将结果与解析解进行比较。 先建立一个该函数的m文件fxy1.m： function f=f(x,y) f=-2*y./x+4*x %注意使用点运算符 return 再输入命令： X,Y=ode45(fxy1,1,3,2); X %显示自变量的一组采样点 Y %显示求解函数与采样点对应的一组数值解 (X.2+1./X.2) %显示求解函数与采样点对应的一组解析解,例6.27 求解初值问题在区间0，2中的解。 建立一个函数文件 fxy2.m： function f=f(x,y) f(2)=-x.*y(2)+x.2-5; f
37、(1)=y(2); f=f; return 在MATLAB命令窗口，输入命令： X,Y=ode45(fxy2,0,2,5,6); X,Y,6.6 非线性方程的数值求解,1单变量非线性方程求解 MATLAB中，提供了求解单变量方程的函数fzero(f,x0,tol)，该函数采用迭代法计算函数f(x)的一个零点，迭代初值为x0，当两次迭代结果小于tol时停止迭代过程。tol的缺省值是eps。注意，在调用函数fzero 之前，要使用m文件建立自己要计算的函数f(x)，只有定义了函数f(x)的m文件后，才能在fzero函数的参数中使用自定义函数名。,例6.28 求f(x)=x-1/x+5 在x0=-5
38、和x0=1作为迭代初值时的零点。 先编制一个函数文件fz.m： function f=f(x) f=x-1/x+5; 然后，在MATLAB命令窗口，输入命令： fzero(fz,-5) %以-5作为迭代初值 Zero found in the interval: -4.8, -5.2. fzero(fz,1),2. 非线性方程组的求解对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为：X=fsolve(fun,X0,option)其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定。
39、最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。,例6.29 求方程组在(1，1，1)附近的解并对结果进行验证。 首先建立方程的函数文件fxyz1.m： function F=F(X) x=X(1);y=X(2);z=X(3); F(1)=sin(x)+y+z2*exp(x); F(2)=x+y*
40、z; F(3)=x*y*z; 在MATLAB命令窗口，输入命令： X=fsolve(fxyz1,1,1,1) %求解X的三个分量x、y、z Y=fxyz1(X) %检验所求结果X是否满足原方程组 norm(Y) %求Y向量的模,例6.30 求圆和直线的两个交点。 建立方程组函数文件fxyz2.m： function F=F(X) x=X(1);y=X(2);z=X(3); F(1)=x2+y2+z2-9; F(2)=3*x+5*y+6*z; F(3)=x-3*y-6*z-1; 在MATLAB命令窗口，输入命令： X1=fsolve(fxyz2,-1,1,-1) %求直线与球面第一个交点 X2=
41、fsolve(fxyz2,1,-1,1) %求直线与球面的第二个交点,6.7 稀疏矩阵,6.7.1 矩阵存储方式 矩阵的完全存储模式完全存储方式是将矩阵的全部元素按列存储。以前讲到的矩阵的存储方式都是按这个方式存储的，此存储方式对稀疏矩阵也适用. 2. 稀疏矩阵的存储方式稀疏存储方式仅存储矩阵所有的非零元素的值及其位置，即行号和列号。在MATLAB中，稀疏存储方式也是按列存储的。 注意!: 在讲稀疏矩阵时，有两个不同的概念，一是指矩阵的0元素较多，该矩阵是一个具有稀疏特征的矩阵，二是指采用稀疏方式存储的矩阵。,6.7.2 稀疏存储方式的产生与转化 1. 将一个完全存储方式的转化为稀疏存储方式
42、函数B=sparse(A)将矩阵A转化为稀疏存储方式的矩阵B。 sparse函数还有其他一些格式： sparse(m,n) 生成一个mn的所有元素都是0的稀疏矩阵。 sparse(u,v,S) u、v、S是三个等长的向量。 此外，还有一些和稀疏矩阵操作有关的函数。例如 U,V,S=find(A) 返回矩阵A中非0元素的下标和元素。这里产生的U、V、S可作为sparse(u,v,s)的参数。 full(A) 返回和稀疏存储矩阵A对应的完全存储方式矩阵。,2. 产生一个稀疏矩阵 把要建立的稀疏矩阵的非0元素及其所在行和列的位置表示出来后由MATLAB自己产生其稀疏存储方式，这需要使用spconve
43、rt函数。调用格式为：B=spconvert(A)其中A为一个m3或m4的矩阵，其每行表示一个非0元素，m是非0元素的个数。 3. 单位稀疏矩阵的产生 单位矩阵只有对角线元素为1，其他元素都为0，是一种具有稀疏特征的矩阵。我们知道，函数eye产生一个完全存储方式的单位矩阵。MATLAB还有一个产生稀疏存储方式的单位矩阵的函数，这就是speye。函数speye(m,n)返回一个mn的稀疏存储单位矩阵。,6.7.3 稀疏矩阵应用举例 稀疏存储矩阵只是矩阵的存储方式不同，它的运算规则与普通矩阵是一样的。所以，在运算过程中，稀疏存储矩阵可以直接参与运算。当参与运算的对象不全是稀疏存储矩阵时，所得结果一般是完全存储形式。,
