《6MATLAB数值计算》由会员分享，可在线阅读，更多相关《6MATLAB数值计算（103页珍藏版）》请在金锄头文库上搜索。
1、第六章 MATLAB数值计算o 数据处理与多项式计算o 数值微积分o 离散傅里叶变换o 线性方程组求解o非线性方程与最优化问题o常微分方程的数值求解o稀疏矩阵2021/7/116.1 数据处理与多项式计算数据处理与多项式计算 6.1.1 数据统计与分析1. 求矩阵最大元素和最小元素MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min，两个函数的调用格式和操作过程类似。（1）求向量的最大值和最小值 y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。 y,I=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取
2、最大值。求向量X的最小值的函数是min(X)，用法和max(X)完全相同。2021/7/12例 求向量x的最大值。命令如下：x=-43,72,9,16,23,47;y=max(x) %求向量x中的最大值y= 72y,l=max(x) %求向量x中的最大值及其该元素的位置y= 72l= 2以上是对行向量进行操作，事实上对列向量的操作与对行向量的操作结果是一样的。例如，对上述x做一转置，有相同的结果：y,l=max(x)y= 722021/7/13（2）求矩阵的最大值和最小值求矩阵A的最大值的函数有3种调用格式，分别是： max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值。
3、Y,U=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。 max(A,dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。 求最小值的函数是min，其用法和max完全相同。2021/7/14例6.1 分别矩阵A中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。A=13,-56,78;25,63,-235;78,25,563;1,0,-1;max(A,2) %求每行最大元素min(A,2) %求每行最小元素max(A) %求每列最大元素min(A) %求每
4、列最小元素max(max(A) %求整个矩阵的最大元素。也可使用命令：max(A(:)min(min(A) %求整个矩阵的最小元素。也可使用命令：min(A(:)2021/7/15（3）两个向量或矩阵对应元素的比较 函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为： 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所有同一位置上的
5、较大元素构成的新矩阵p。2021/7/16 x=4,5,6;1,4,8x = 4 5 6 1 4 8 y=1,7,5;4,5,7y = 1 7 5 4 5 7 p=max(x,y) %在x,y同一位置上的两个元素中找出较大值p = 4 7 6 4 5 8以上是对两个大小的矩阵操作，MATLAB还允许对一个矩阵和一个常数或单变量操作。例如，仍然用上例的矩阵x和已赋值为4.5的变量f，操作如下： f=4.5; p=max(x,f)p = 4.5000 5.0000 6.0000 4.5000 4.5000 8.00002021/7/172. 求矩阵的平均值和中值求数据序列平均值的函数是mean，求
6、数据序列中值的函数是median。两个函数的调用格式为：mean(X)：返回向量X的算术平均值。median(X)：返回向量X的中值。mean(A)：返回一个行向量，其第i个元素是A的第i列的算术平均值。median(A)：返回一个行向量，其第i个元素是A的第i列的中值。mean(A,dim)：当dim为1时，该函数等同于mean(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的算术平均值。median(A,dim)：当dim为1时，该函数等同于median(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的中值。2021/7/18例：求向量x的平均值和中值，操作
7、如下： y=9,-2,5,6,7,12; mean(y)ans = 6.1667 median(y)ans = 6.50002021/7/193. 矩阵元素求和与求积数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，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个元素是A的第i
8、行的各元素之和。prod(A,dim)：当dim为1时，该函数等同于prod(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。2021/7/110例6.2 求矩阵A的每行元素的乘积和全部元素的乘积。A=1,2,3,4;5,6,7,8;9,10,11,12;S=prod(A,2)S = 24 1680 11880prod(S) %求A的全部元素的乘积。也可以使用命令prod(A(:) ans = 4790016002021/7/1114. 矩阵元素累加和与累乘积在MATLAB中，使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的
9、调用格式为：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,dim)：当dim为1时，该函数等同于cumprod(A)；当dim为2时，返回一个向量，其第i行是A的第i行的累乘积向量。2021/7/112例例6.3 求向量X=(1！，2！，3！，10！)。X=cu
10、mprod(1:10)X = Columns 1 through 6 1 2 6 24 120 720 Columns 7 through 10 5040 40320 362880 36288002021/7/113对于具有N个元素的数据序列 ，标准方差的计算公式如下：2021/7/1145求标准方差在MATLAB中，提供了计算数据序列的标准方差的函数std。对于向量X，std(X)返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。std函数的一般调用格式为：Y=std(A,flag,dim)其中dim取1或2。当dim=1时，求各列元素的标
11、准方差；当dim=2时，则求各行元素的标准方差。flag取0或1，当flag=0时，按S1所列公式计算标准方差，当flag=1时，按S2所列公式计算标准方差。缺省flag=0，dim=1。2021/7/115例6.4 对二维矩阵x，从不同维方向求出其标准方差。x=4,5,6;1,4,8 % 产生一个二维矩阵xy1=std(x,0,1)y2=std(x,1,1)y3=std(x,0,2)y4=std(x,1,2)2021/7/1166相关系数MATLAB提供了corrcoef函数，可以求出数据的相关系数矩阵。corrcoef函数的调用格式为：corrcoef(X)：返回从矩阵X形成的一个相关系数
12、矩阵。此相关系数矩阵的大小与矩阵X一样。它把矩阵X的每列作为一个变量，然后求它们的相关系数。corrcoef(X,Y)：在这里，X,Y是向量，它们与corrcoef(X,Y)的作用一样。2021/7/117例6.5 生成满足正态分布的100005随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)2021/7/1187. 排序MATLAB中对向量X是排序函数是sort(X)，函数返回一个对X中的元素按升序排列的新向量。sort函数也可以对矩阵A的各列或各行重新排序，其
13、调用格式为：Y,I=sort(A,dim)其中dim指明对A的列还是行进行排序。若dim=1，则按列排；若dim=2，则按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。2021/7/119例例6.6 对下列矩阵做各种排序。 A=1,-8,5;4,12,6;13,7,-13;sort(A) %对A的每列按升序排序-sort(-A,2) %对A的每行按降序排序 X,I=sort(A) %对A按列排序，并将每个元素所在行号送矩阵I2021/7/1206.1.2 数据插值在工程测量和科学实验中，所得到的数据通常都是离散的。如果要得到这些离散点以外的其他点的数值，就需要根据这些已知数据进行插值。
14、例如，测量得n个点的数据 ，这些数据点反映了一个函数关系 ，然而并不知道 f(x) 的解析式。数据插值的任务就是根据上述条件构造一个函数 使得对于 有 ，且在两个相邻采样点 ，g(x)光滑过渡。插值函数一般由线性函数、多项式、样条函数或这些函数的分段函数充当。2021/7/1216.1.2 数据插值1. 一维数据插值在MATLAB中，实现这些插值的函数是interp1，其调用格式为：Y1=interp1(X,Y,X1,method)函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。
15、method是插值方法，允许的取值有linear、nearest、cubic、spline。2021/7/122注意：X1的取值范围不能超出X的给定范围，否则，会给出“NaN”错误。 例6.7 给出概率积分的数据表如表6.1所示，用不同的插值方法计算f(0.472)。 表6.1 概率积分数据表x0.460.470.480.49f(x)0.48465550.49375420.50274980.51166832021/7/123x=0.46:0.01:0.49; %给出x，f(x)f=0.4846555,0.4937542,0.5027498,0.5116683;format longinterp
16、1(x,f,0.472) %用默认方法，即线性插值方法计算f(x)interp1(x,f,0.472,nearest) %用最近点插值方法计算f(x)interp1(x,f,0.472,spline) %用3次样条插值方法计算f(x)interp1(x,f,0.472,cubic) %用3次多项式插值方法计算f(x)format short2021/7/124例例6.8 某检测参数f随时间t的采样结果如表6.2，用数据插值法计算t=2，7，12，17，22，17，32，37，42，47，52，57时的f值。表6.2 检测参数f随时间t的采样结果t051015202530f3.10252.256
17、879.51835.92968.84136.25237.9t35404550556065f6152.76725.36848.36403.56824.77328.57857.62021/7/125T=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,7328.5,7857.6;F1=interp1(T,F,X) %用线性插值方法插值F1=interp1(T,F,X,nearest) %用最近点插值方法插值F1=interp1(T,F,X,spl
18、ine) %用3次样条插值方法插值F1=interp1(T,F,X,cubic) %用3次多项式插值方法插值2021/7/1262. 二维数据插值在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式为：Z1=interp2(X,Y,Z,X1,Y1,method)其中X,Y是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的函数值，X1,Y1是两个向量或标量，描述欲插值的点。Z1是根据相应的插值方法得到的插值结果。 method的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。同样，X1,Y1的取值范围不能超出X,Y的给定范围，否则，会给出“NaN”错误。2021
19、/7/127例例6.9 设z=x2+y2，对z函数在0,10,2区域内进行插值。x=0:0.1:1;y=0:0.2:2;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)点插值%下一命令在(0.5,0.4),(0.6,0.4),(0.5,0.5)和(0.6,0.5)各点插
20、值interp2(x,y,Z,0.5 0.6,0.4 0.5)2021/7/128例例6.10 某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点(米)，用h表示测量时间(秒)，用T表示测得各点的温度()，测量结果如表6.2所示表6.2 钢轨各点温度测量值 T xh0 2.5 5 7.5 10095 14 0 0 03088 48 32 12 66067 64 54 48 41试用用3次多项式插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。2021/7/129x=0:2.5:10;h=0:30:60;T=95,14,0,0,0;88,48,32,12,6;67,64,54
21、,48,41;xi=0:0.5:10;hi=0:10:60;temps=interp2(x,h,T,xi,hi,cubic);mesh(xi,hi,temps);2021/7/1306.1.3 曲线拟合与数值插值类似，曲线拟合的目的也是用一个较简单的函数去逼近一个复杂的或者未知的函数，所依据的条件都是在一个区间或一个区域上的有限个采样点的函数值。数值插值要求逼近函数在采样点与被逼近函数相等，但由于实验或测量中的误差，所获得的数据不一定准确。在这种情况下，如果强求逼近函数通过采样点，显然是不够合理的。为此，构造函数y=g(x)去逼近f(x) ，这里不要求曲线g(x) 严格通过采样点，但希望g(x
22、)能尽量地靠近这些点，就是使误差g(xi)-f(xi) 在某种意义上达到最小。MATLAB曲线拟合的最优标准是采用常见的最小二乘原理，所构造的g(x)是一个次数小于插值节点个数的多项式。2021/7/1316.1.3 曲线拟合在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。polyfit函数的调用格式为：P,S=polyfit(X,Y,m)函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数。polyva
23、l函数的功能是按多项式的系数计算x点多项式的值。2021/7/132例例6.11 用一个3次多项式在区间0，2内逼近函数 X=linspace(0,2*pi,50);Y=sin(X);P=polyfit(X,Y,3) %得到3次多项式的系数和误差以上求得了3次拟合多项式p(x)的系数，得p(x)=0.0912x3-0.8596x2+1.8527x-0.1649X=linspace(0,2*pi,20);Y=sin(X);Y1=polyval(P,X)plot(X,Y,:o,X,Y1,-*)2021/7/1336.1.4 多项式计算在MATLAB中，n次多项式用一个长度为n+1的行向量表示，缺少
24、的幂次项系数为0.如：则在MATLAB中，p(x)表示为向量形式：1. 多项式的四则运算（1）多项式的加减运算（2）多项式乘法运算函数conv(P1,P2)用于求多项式P1和P2的乘积。这里，P1、P2是两个多项式系数向量。2021/7/134（3）多项式除法函数Q,r=deconv(P1,P2)用于对多项式P1和P2作除法运算。其中Q返回多项式P1除以P2的商式，r返回P1除以P2的余式。这里，Q和r仍是多项式系数向量。deconv是conv的逆函数，即有P1=conv(P2,Q)+r。例例6.12 设(1)求f(x)+g(x)、f(x)-g(x)。(2)求f(x)g(x)、f(x)/g(x
25、)。2021/7/135f=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) %求f(x)/g(x)，商式送Q，余式送r。 2021/7/1362. 多项式的导函数对多项式求导数的函数是：p=polyder(P)：求多项式P的导函数p=polyder(P,Q)：求PQ的导函数p,q=polyder(P,Q)：求P/Q的导函数，导函数的分子存入p，分母存入q。上述函数中，参数P,Q是多项式的向量表示，结果p,q也是多项式的向量表示。2021
26、/7/137例例6.13 求有理分式的导数。P=3,5,0,-8,1,-5;Q=10,5,0,0,6,0,0,7,-1,0,-100;p,q=polyder(P,Q)2021/7/1383. 多项式求值MATLAB提供了两种求多项式值的函数：polyval与polyvalm，它们的输入参数均为多项式系数向量P和自变量x。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。(1)代数多项式求值polyval函数用来求代数多项式的值，其调用格式为：Y=polyval(P,x)若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。2021/7/139
27、例例6.14 已知多项式x4+8x3-10，分别取x=1.2和一个23矩阵为自变量计算该多项式的值。A=1,8,0,0,-10; % 4次多项式系数x=1.2; % 取自变量为一数值y1=polyval(A,x)x=-1,1.2,-1.4;2,-1.8,1.6 % 给出一个矩阵xy2=polyval(A,x) % 分别计算矩阵x中各元素为自变量的多项式之值2021/7/140(2)矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。设A为方阵，P代表多项式x3-5x2+8，那么pol
28、yvalm(P,A)的含义是：A*A*A-5*A*A+8*eye(size(A)而polyval(P,A)的含义是：A.*A.*A-5*A.*A+8*ones(size(A)例6.15 仍以多项式x4+8x3-10为例，取一个22矩阵为自变量分别用polyval和polyvalm计算该多项式的值。2021/7/141A=1,8,0,0,-10; % 多项式系数x=-1,1.2;2,-1.8 % 给出一个矩阵xy1=polyval(A,x) % 计算代数多项式的值y2=polyvalm(A,x) % 计算矩阵多项式的值4. 多项式求根n次多项式具有n个根，当然这些根可能是实根，也可能含有若干对共
29、轭复根。MATLAB提供的roots函数用于求多项式的全部根，其调用格式为：x=roots(P)其中P为多项式的系数向量，求得的根赋给向量x，即x(1),x(2),x(n)分别代表多项式的n个根。2021/7/142例6.16 求多项式x4+8x3-10的根。命令如下：A=1,8,0,0,-10;x=roots(A)若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：P=poly(x)若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。2021/7/143例6.17 已知 f(x)(1) 计算f(x)=0 的全部根。(2) 由方程
30、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)2021/7/1446.2 数值微积分数值微积分6.2.1 数值微分一般来说，函数的导数依然是一个函数。设函数f(x)的导函数f(x)=g(x),高等数学关心的是g(x)的形式和性质，而数值分析关心的问题是怎样计算g(x)在一串离散点的近似值以及所计算的近似值有多大误差。1.数值差分与差商任意函数f(x)在x点的导数是通过极限定义的2021/7/145上述式子中，均假设 ，如果去掉上述等式右端的 的极限过
31、程，并引进记号：称 及 分别为函数在x点处以h为步长的向前差分、向后差分和中心差分。当步长h充分小时，有和差分一样，称 及 分别为函数在以h为步长的向前差商、向后差商和中心差商。2021/7/146当步长h充分小时，函数f在x点的微分接近于函数在该点的任意种差分，而f在x点的导数接近于函数在该点的任意种差商。2. 数值微分的实现在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数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)=d
32、iff(diff(X)。DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。2021/7/147例6.18 设x由0,2间均匀分布的10个点组成，求sinx的13阶差分。命令如下：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)2021/7/148例6.19 用不同的方法求函数f(x)的数值导数，
33、并在同一个坐标系中做出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.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的导函数
34、g在假设点的导数plot(x,dpx,x,dx,.,x,gx,-); %作图2021/7/1496.2.2 数值积分1. 数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。基本梯形与辛普林求积公式复合梯形与辛普林求积公式2021/7/1502. 数值积分的实现被积函数一般是用一个解析式给出，但也有很多情况下用一个表格给出。在MATLAB中，对这两种
35、给定被积函数的方法，提供了不同的数值积分函数。(1)被积函数是一个解析式MATLAB提供了quad函数和quadl函数来求定积分。它们的调用格式为： quad(filename,a,b,tol,trace) quadl(filename,a,b,tol,trace)其中filename是被积函数名。a,b分别是定积分的下限和上限。tol用来控制积分精度，默认取10-6.trace控制是否展现积分过程，非0展现，取0不展现，默认取0.2021/7/151例6.20 用两种不同的方法求定积分。先建立一个函数文件ex.m：function ex=ex(x)ex=exp(-x.2);然后在MATLAB
36、命令窗口，输入命令：format longI=quad(ex,0,1) %注意函数名应加字符引号I = 0.74682418072642I=quadl(ex,0,1)I = 0.74682413398845也可不建立关于被积函数的函数文件，而使用语句函数(内联函数)求解，命令如下：g=inline(exp(-x.2); %定义一个语句函数g(x)=exp(-x2)I=quadl(g,0,1) %注意函数名不加号I = 0.74682413398845format short2021/7/152(2)被积函数由一个表格定义在科学实验和工程应用中，函数关系往往是不知道的，只有实验测定的一组样本点和
37、样本值，这时，就无法使用quad函数计算其定积分。在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X、Y定义函数关系Y=f(X)。X、Y是两个等长的向量：X=(x1,x2,xn)，Y=(y1,y2,yn)，并且x1x2 A=1,-1,1;5,-4,3;2,1,1A = 1 -1 1 5 -4 3 2 1 1 L,U=lu(A) LU=L*U利用第二种格式对矩阵A进行LU分解： L,U,P=lu(A) LU=L*U inv(P)*L*P2021/7/166例6.25 用LU分解求解例6.24中的线性方程组。命令如下：A=2,1,-5,1;1,-5,0
38、,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)2021/7/167 (2) QR分解对矩阵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分解后，线性方程组A
39、x=b的解x=R(Qb)或x=E(R(Qb)。2021/7/168设：则对矩阵A进行QR分解的命令如下：A1,-1,1;5,-4,3;2,7,10；Q,R=qr(A)为检验结果是否正确，输入命令：QR=Q*R利用第二种格式对矩阵A进行QR分解Q,R,E=qr(A)Q*R/E%验证A= Q*R*inv(E)2021/7/169例6.26 用QR分解求解例6.24中的线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;Q,R=qr(A);x=R(Qb)或采用QR分解的第2种格式，命令如下：Q,R,E=qr(A);x=E*(R(
40、Qb)将得到与上面同样的结果2021/7/170 (3) Cholesky分解如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即X=RR。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为：R=chol(X)：产生一个上三角阵R，使RR=X。若X为非对称正定，则输出一个出错信息。R,p=chol(X)：这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足RR=X(1
41、:q,1:q)。实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x=R(Rb)。2021/7/171设则对矩阵A进行cholesky分解的命令如下：A=2,1,1;1,2,-1;1,-1,3;R=chol(A)可以验证RR=A:R*R利用第二种格式对矩阵A进行cholesky分解：R,P=chol(A)p=0表示A是一个正定阵，对一个非正定矩阵会给出错误信息，所以，chol函数可用来判定矩阵是否为正定矩阵2021/7/172例6.27 用Cholesky分解求解例6.24中的线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b
42、=13,-9,6,0;R=chol(A)? Error using = cholMatrix must be positive definite命令执行时，出现错误信息，说明A为非正定矩阵。2021/7/1736.4.2 迭代解法迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。为了求解线性方程组好处是将一组x代入右端，可以立即得到另一组x,若两组x相等，那么它就是方程组的解，不等时可以继续迭代。2021/7/1741Jacobi迭代法对于线性方程组Ax=b，如果A为非奇异方阵，即ai
43、i0(i=1,2,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是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的解。2021/7/175Jacobi迭代法的迭代法的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
44、+f; n=n+1;end2021/7/176例6.28 用Jacobi迭代法求解线性方程组。设迭代初值为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)2021/7/1772Gauss-Serdel迭代法在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：x(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为Ga
45、uss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。2021/7/178Gauss-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;end2021/7/179例6.29 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在命令中调用函数文件gauseidel.m，
46、命令如下：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)注：一般情况下Gauss迭代比Jacobi快。但也不是绝对的，某些情况下，Jacobi收敛而Gauss却可能不收敛，看下例：2021/7/180例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=gauseidel(a,b,0;0;0)2021/7/1816.4.3 求线性方程组的通解求线性方
47、程组的通解线性方程组的求解分为两类：一类是求方程组的惟一解即特解，另一类是求方程组的无穷解即通解。这里对线性方程组 Ax=b的求解理论作一个归纳。(1)当系数矩阵A是一个满秩方阵时，方程Ax=b称为恰定方程，方程有惟一解x=A-1b，这是最基本的一种情况。一般用x=Ab求解速度更快。(2)当方程组右端向量b=0时，方程称为齐次方程组。齐次方程组总有零解，因此称解x=0为平凡解。当系数矩阵A的秩小于n(n为方程组中未知变量的个数)时，齐次方程组有无穷多个非平凡解，其通解中包含n-rank(A)个线性无关的解向量，用MATLAB的函数null(A,r)可求得基础解系。2021/7/182(3)当方
48、程组右端向量b0时，系数矩阵的秩rank(A)与其增广矩阵的秩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(A)个线性无关的解向量。当rank(A)0 %非齐次方程组非齐次方程组if rank(A)=rank(A,b)if rank(A)=n %有惟一解有惟一
49、解disp(原方程组有惟一解原方程组有惟一解x);x=Ab;else %方程组有无穷多个解，基础解系方程组有无穷多个解，基础解系disp(原方程组有无穷个解，特解为原方程组有无穷个解，特解为x，其齐次方程组的基础解系为，其齐次方程组的基础解系为y);x=Ab; y=null(A,r);endelsedisp(方程组无解方程组无解); %方程组无解方程组无解x=;endelse %齐次方程组齐次方程组disp(原方程组有零解原方程组有零解x);x=zeros(n,1); %0解解if rank(A)n disp(方程组有无穷个解，基础解系为方程组有无穷个解，基础解系为y); %非非0解解y=nu
50、ll(A,r);endend2021/7/184例例6.31 求解方程组A=1,-2,3,-1;3,-1,5,-3;2,1,2,-2;b=1;2;3;x,y=line_solution(A,b)2021/7/185例例6.32 求方程组的通解format rat %指定有理式格式输出A=1,1,-3,-1;3,-1,-3,4;1,5,-9,-8;b=1,4,0;x,y=line_solution(A,b);x,yformat short %恢复默认的短格式输出说明：原方程组有无穷多个解，特解为x，其齐次方程组的基础解系为y2021/7/1866.5 非线性方程与最优化问题求解非线性方程与最优化
51、问题求解6.5.1 非线性方程数值求解非线性方程的求根方法很多，常用的有牛顿迭代法，但该方法需要求原方程的导数，而在实际运算时这一条件有时是不能满足的，所以又出现了弦截法、二分法等其他方法。1. 单变量非线性方程求解 在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。该函数的调用格式为： z=fzero(fname,x0,tol,trace)其中fname是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示
52、，缺省时取trace=0。2021/7/187例6.33 求 在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 2021/7/1882. 非线性方程组的求解 对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为： X=fsolve(fun,X0,option)其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件
53、名，X0是求根过程的初值，option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。2021/7/189例6.34 求下列方程组在(1，1，1)附近的解并对结果进行验证。首先建立函数文件myfun.m。function F=myfun (X)x=X(1);y=X(2);
54、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下，调用fsolve函数求方程的根。X=fsolve(myfun,1,1,1,optimset(Display, off)X = 0.0224 -0.0224 -0.0000 将得到的解代回原方程，可以检验结果是否正确，命令为：q=myfun(X)2021/7/190例例6.35 求圆和直线的两个交点。圆：直线：先建立方程组函数文件fxyz.m：function F=fxyz(X)x=X(1);y=X(2);z=X(3);F(1)=x2+y2+z2-
55、9;F(2)=3*x+5*y+6*z;F(3)=x-3*y-6*z-1;再在MATLAB命令窗口，输入命令：X1=fsolve(fxyz,-1,1,-1,optimset(Display, off) %求第一个交点X2=fsolve(fxyz,1,-1,1,optimset(Display, off) %求第二个交点使用使用fsolve函数求解方程组函数求解方程组时，必须先估计出方程组时，必须先估计出方程组的根的大致范围。的根的大致范围。2021/7/1916.5.2 无约束最优化问题求解在实际应用中，许多科学研究和工程计算问题都可以归结为一个最小化问题，如能量最小、时间最短等。MATLAB提
56、供了3个求最小值的函数，它们的调用格式为：(1)x,fval=fminbnd(filename,x1,x2,option)：求一元函数在(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)上的最小值就是f(x)在(a,b)的最大值，所以fminbnd(-f
57、,x1,x2)返回函数f(x)在区间(x1,x2)上的最大值。2021/7/192例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)内的最小值点和最小值fminbnd(f,1,10) %求函数在(1，10)内的最小值点。注意函数名f不用加例6.37 求函数 在(0.5,0.5,0.5)附近的最小值。建立函数文件fxyz.m：func
58、tion f=fxyz(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;在MALAB命令窗口，输入命令：U,fmin=fminsearch(fxyz,0.5,0.5,0.5) %求函数的最小值点和最小值2021/7/1936.5.3 有约束最优化问题求解约束条件可以进一步细化为：线性不等式约束： 线性等式约束：非线性不等式约束： 非线性等式约束：x的下界和上界：MATLAB最优化工具箱提供了一个fmincon函数，专门用于求解各种约束下的最优化问题。该函数的调用格式为：x,fval=fmincon(filename,x0,A,b, Aeq,beq
59、,Lbnd,Ubnd, NonF,option)其中x、fval、filename、x0和option的含义与求最小值函数相同。其余参数为约束条件，参数NonF为非线性约束函数的M文件名。如果某个约束不存在，则用空矩阵来表示。2021/7/194例6.38 求解有约束最优化问题。首先编写目标函数M文件fop.m。function f=fop(x)f=0.4*x(2)+x(1)2+x(2)2-x(1)*x(2)+1/30*x(1)3;再设定约束条件，并调用fmincon函数求解此约束最优化问题。x0=0.5;0.5;A=-1,-0.5;-0.5,-1;b=-0.4;-0.5;lb=0;0;opt
60、ion=optimset; option.LargeScale=off; option.Display =off;x,f=fmincon(fop ,x0,A,b,lb,option)2021/7/1956.6 常微分方程的数值求解常微分方程的数值求解只有对一些典型的常微分方程，才能求出它们的一般解表达式并用初始条件确定表达式中的任意常数。然而在实际问题中遇到的常微分方程往往很复杂，在许多情况下得不出一般解，所以，一般是要求获得解在若干个点上的近似值。考虑常微分方程的初值问题：所谓其数值解法，就是求它的解y(t)在节点处的近似值的方法。所求得的近似值称为常微分方程初值问题的数值解。一般采用等距节
61、点，h是步长。常用的有欧拉方法，龙格-库塔法，线性多步法，预报校正法2021/7/1966.6 常微分方程的数值求解常微分方程的数值求解6.6.1 龙格库塔法简介龙格库塔法简介龙格库塔公式由中值定理而来。其中：2021/7/1976.6.2 龙格库塔法的实现 基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0)其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为t0,tf,表示求解区间。y0是初始状态列向量。t和y分别给出时间向
62、量和相应的状态向量。这两个函数分别采用了二阶、三阶龙格-库塔法和四阶、五阶龙格-库塔法，并采用自适应变步长的求解方法，即当解的变化较慢时采用较大的步长，从而使得计算速度很快，当解的变化较快时步长会自动地变小，从而使得计算精度很高。2021/7/198例6.39 设有初值问题，试求其数值解，并与精确解相比较(精确解为y(t)= )。 (1) 建立函数文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。t0=0;tf=10;y0=2;t,y=ode23(funt,t0,tf,y0); %求数值解y1=sqrt(t+1)+1; %
63、求精确解plot(t,y,b.,t,y1,r-)y(0)=22021/7/199例例6.40 已知一个二阶线性系统的微分方程为：其中a=2，绘制系统的时间响应曲线和相平面图。建立一个函数文件sys.m：function xdot=sys(t,x)xdot=-2*x(2);x(1); 取t0=0，tf=20，求微分方程的解：t0=0;tf=20;t,x=ode45(sys,t0,tf,1,0); t,xsubplot(1,2,1);plot(t,x(:,2); %解的曲线，即t-xsubplot(1,2,2);plot(x(:,2),x(:,1) %相平面曲线，即x-xaxis equal2021/7/11006.7 稀疏矩阵稀疏矩阵6.7.1 矩阵存储方式矩阵存储方式MATLAB的矩阵有两种存储方式：完全存储方式和稀疏存储方式。1完全存储方式完全存储方式是将矩阵的全部元素按列存储。以前讲到的矩阵的存储方式都是按这个方式存储的，此存储方式对稀疏矩阵也适用。 2稀疏存储方式稀疏存储方式仅存储矩阵所有的非零元素的值及其位置，即行号和列号。在MATLAB中，稀疏存储方式也是按列存储的。2021/7/11012021/7/1102 结束语结束语若有不当之处，请指正，谢谢！若有不当之处，请指正，谢谢！
