第7章 MATLAB在计算方法中应用

1、第第7章章 MATLAB在计算方法中的应用在计算方法中的应用6.1 数据统计处理数据统计处理6.1.1 最大值和最小值最大值和最小值MATLAB提供的求数据序列的最大值和最小提供的求数据序列的最大值和最小值的函数分别为值的函数分别为max和和min，两个函数的调，两个函数的调用格式和操作过程类似。用格式和操作过程类似。1求向量的最大值和最小值求向量的最大值和最小值求一个向量求一个向量X的最大值的函数有两种调用格的最大值的函数有两种调用格式，分别是：式，分别是：(1) y=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，如果如果X中包含复数元素，则按模取最大值。中包含复数元素，则按

2、模取最大值。(2) y,I=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，最大，最大值的序号存入值的序号存入I，如果，如果X中包含复数元素，则按模中包含复数元素，则按模取最大值。取最大值。求向量求向量X的最小值的函数是的最小值的函数是min(X)，用法和，用法和max(X)完全相同。完全相同。例例6-1 求向量求向量x的最大值。的最大值。命令如下：命令如下：x=-43,72,9,16,23,47;y=max(x) %求向量求向量x中的最大值中的最大值y,l=max(x) %求向量求向量x中的最大值及其该元素中的最大值及其该元素的位置的位置2求矩阵的最大值和最小值求矩阵的最大值和

3、最小值求矩阵求矩阵A的最大值的函数有的最大值的函数有3种调用格式，分种调用格式，分别是：别是：(1) max(A)：返回一个行向量，向量的第：返回一个行向量，向量的第i个元个元素是矩阵素是矩阵A的第的第i列上的最大值。列上的最大值。(2) Y,U=max(A)：返回行向量：返回行向量Y和和U，Y向量向量记录记录A的每列的最大值，的每列的最大值，U向量记录每列最向量记录每列最大值的行号。大值的行号。(3) max(A,dim)：dim取取1或或2。dim取取1时，时，该函数和该函数和max(A)完全相同；完全相同；dim取取2时，该时，该函数返回一个列向量，其第函数返回一个列向量，其第i个元素是

4、个元素是A矩阵矩阵的第的第i行上的最大值。行上的最大值。求最小值的函数是求最小值的函数是min，其用法和，其用法和max完全相完全相同。同。例例6-2 分别求矩阵分别求矩阵x中各列和各行元素中的最中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。大值，并求整个矩阵的最大值和最小值。 A=magic(4)3两个向量或矩阵对应元素的比较两个向量或矩阵对应元素的比较函数函数max和和min还能对两个同型的向量或矩阵还能对两个同型的向量或矩阵进行比较，调用格式为：进行比较，调用格式为：(1) U=max(A,B)：A,B是两个同型的向量或矩是两个同型的向量或矩阵，结果阵，结果U是与是与A,B同

5、型的向量或矩阵，同型的向量或矩阵，U的的每个元素等于每个元素等于A,B对应元素的较大者。对应元素的较大者。(2) U=max(A,n)：n是一个标量，结果是一个标量，结果U是与是与A同型的向量或矩阵，同型的向量或矩阵，U的每个元素等于的每个元素等于A对应对应元素和元素和n中的较大者。中的较大者。min函数的用法和函数的用法和max完全相同。完全相同。例例6-3 求两个矩阵求两个矩阵x, y所有同一位置上的较大元所有同一位置上的较大元素构成的新矩阵素构成的新矩阵p。B=pascal(4)6.1.2 求和与求积求和与求积数据序列求和与求积的函数是数据序列求和与求积的函数是sum和和prod，其使用

6、方法类似。设其使用方法类似。设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时，返回一个列向量，时，返回一

7、个列向量，其第其第i个元素是个元素是A的第的第i行的各元素之和。行的各元素之和。prod(A,dim)：当：当dim为为1时，该函数等同于时，该函数等同于prod(A)；当；当dim为为2时，返回一个列向量，时，返回一个列向量，其第其第i个元素是个元素是A的第的第i行的各元素乘积。行的各元素乘积。例例6-4 求矩阵求矩阵A的每行元素的乘积和全部元素的每行元素的乘积和全部元素的乘积。的乘积。6.1.3 平均值和中值平均值和中值求数据序列平均值的函数是求数据序列平均值的函数是mean，求数据序列中值的函数是，求数据序列中值的函数是median。两个函数的调用格式为：。两个函数的调用格式为：mean

8、(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)：当：当di

9、m为为1时，该函数等同于时，该函数等同于median(A)；当；当dim为为2时，返回一个列向量，其第时，返回一个列向量，其第i个元素是个元素是A的第的第i行的中行的中值。值。例例6-5 求向量求向量x的平均值和中值。的平均值和中值。x=A(:,1) mean(x) median(x)6.1.4 累加和与累乘积累加和与累乘积在在MATLAB中，使用中，使用cumsum和和cumprod函数能方便地求得函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式向量和矩阵元素的累加和与累乘积向量，函数的调用格式为：为：cumsum(X)：返回向量：返回向量X累加和向量。累加和向量。cump

10、rod(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)；当

11、当dim为为2时，返回一个向量，其第时，返回一个向量，其第i行是行是A的第的第i行的累乘行的累乘积向量。积向量。例例6-6 求向量求向量x累加与累乘积。累加与累乘积。6.1.6 排序排序MATLAB中对向量中对向量X是排序函数是是排序函数是sort(X)，函数返，函数返回一个对回一个对X中的元素按升序排列的新向量。中的元素按升序排列的新向量。sort函数也可以对矩阵函数也可以对矩阵A的各列或各行重新排序，其的各列或各行重新排序，其调用格式为：调用格式为：Y,I=sort(A,dim)其中其中dim指明对指明对A的列还是行进行排序。若的列还是行进行排序。若dim=1，则按列排；若则按列排；若di

12、m=2，则按行排。，则按行排。Y是排序后的矩是排序后的矩阵，而阵，而I记录记录Y中的元素在中的元素在A中位置。中位置。sort(x)sort(A)Y,I=sort(A) 7.1 插值插值7.1.1 Lagrange 插值插值 对给定的对给定的n个自变量个自变量x1, x2, , xn及对应的及对应的 函数值函数值y1, y2, , yn, 求插值区间内任意求插值区间内任意x的函的函数值数值y。 lagrange 插值公式：插值公式： y(x)= yk（ ）nk=1nj=1jkx-xjxk-xj将上述插值公式编写成将上述插值公式编写成M文件，取名文件，取名lagrange.m句法结构：句法结构：

13、function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end ends=p*y0(k)+s;endy(i)=s;end例例1：x=0:0.2:2*pi; y=sin(x);x1=0:0.1:2*pi; 插值：插值： y1=lagrange(x,y,x1); 精确值：精确值： y2=sin(x1); 检查插值效果：检查插值效果：plot(x,y,r,x1,y1,bo) plot(x

14、1,y1-sin(x1)例例2：x=-5:1:5; y=1./(1+x.2);x0=-5:0.1:5; 插值：插值： y0=lagrange(x,y,x0); 精确值：精确值： y1= 1./(1+x0.2); 检查插值效果：检查插值效果：plot(x0,y0,-r,x0,y1,-b) 7.1.2 一维插值命令一维插值命令 interp1 句法结构：句法结构： Y1=interp1(X,Y,X1,method) 对一组节点（对一组节点（x,y）进行插值，计算插值点）进行插值，计算插值点 x1的函数值。的函数值。X为节点向量值，为节点向量值，Y为对应的节点为对应的节点函数值。函数值。Y1是一个与

15、是一个与X1等长的插值结果。等长的插值结果。method是插值方法，允许的取值有：是插值方法，允许的取值有：nearest:线性最近项插值，线性最近项插值，linear:线性插值线性插值 spline:三次样条插值，三次样条插值，cubic：三次插值。三次插值。 注意：注意：X1的取值范围不能超出的取值范围不能超出X的给定范的给定范围，否则，会给出围，否则，会给出“NaN”错误。错误。例例 正弦曲线的插值示例。正弦曲线的插值示例。 x=0:0.1:10; y=sin(x); xi=0:0.25:10; yi=interp1(x,y,xi); plot(x,y,o,xi,yi)例例 计算上例的函

16、数插值计算上例的函数插值x=-5:1:5; y=1./(1+x.2);y2=interp1(x,y,x0);plot(x0,y0,-r,x0,y1,-b)hold on;plot(x0,y2,*m)hold off7.1.3 二维数据插值二维数据插值在在MATLAB中，提供了解决二维插值问题的函中，提供了解决二维插值问题的函数数interp2，其调用格式为：，其调用格式为：Z1=interp2(X,Y,Z,X1,Y1,method) 其中其中X,Y是两个向量，分别描述两个参数的采是两个向量，分别描述两个参数的采样点，样点，Z是与参数采样点对应的函数值，是与参数采样点对应的函数值，X1,Y1是两

17、个向量或标量，描述欲插值的点。是两个向量或标量，描述欲插值的点。Z1是根据相应的插值方法得到插的值结果。是根据相应的插值方法得到插的值结果。 method的取值与一维插值函数相同。的取值与一维插值函数相同。X,Y,Z也也可以是矩阵形式。可以是矩阵形式。注意：注意：X1,Y1的取值范围不能超出的取值范围不能超出X,Y的给定范的给定范围，否则，会给出围，否则，会给出“NaN”错误。错误。例：例： 某实验对一根长某实验对一根长10米的钢轨进行热源的温度传米的钢轨进行热源的温度传播测试。用播测试。用x表示测量点表示测量点0:2.5:10(米米)，用，用h表示测表示测量时间量时间0:30:60(秒秒)，

18、用，用T表示测试所得各点的温度表示测试所得各点的温度()。试用线性插值求出在一分钟内每隔。试用线性插值求出在一分钟内每隔20秒、秒、钢轨每隔钢轨每隔1米处的温度米处的温度TI。命令如下：命令如下：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:10;hi=0:20:60;TI=interp2(x,h,T,xi,hi);7.1.4 Hermite插值插值1.方法介绍方法介绍已知已知n个插值节点个插值节点x1,x2,，xn及其对应的函数值及其对应的函数值 y1,y2,yn和一阶导数值和一阶导数值y1,y2,y

19、n。则计算插则计算插值区域内任意值区域内任意x的函数值的函数值y的的Hermite插值公式：插值公式：nijijiinijjjijiiiiiiniixxxxxxhyyyxxhxy12111;)()2)()(其中2 MATLAB实现实现 hermite.m function y=hermite(x0,y0,y1,x) n=length(x0);m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a; end e

20、nd yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i); end y(k)=yy; end例例 Hermite插值示例。对给定数据表如表插值示例。对给定数据表如表7.2所所示，试构造示，试构造Hermite多项式，并给出多项式，并给出sin0.34的近的近似值。似值。 x0=0.3 0.32 0.35; y0=0.29552 0.31457 0.34290; y1=0.95534 0.94924 0.93937; x=0.3:0.005:0.35; y=hermite(x0,y0,y1,x); plot(x,y,x0,y0)y= hermite (x0,y0

21、,y1,0.34); y sin(0.34) 7.1.5 三次样条插值三次样条插值1 方法介绍方法介绍设区间设区间a,b 上给定有关划分上给定有关划分 a=x0 x1=x0(j)&(x(k) x=28.7 28 29 30; y=4.1 4.3 4.1 3.0; x0=28.7:0.15:30; y=spline2(x,y,0,0,x0); plot(x0,y)7.1.6 最小二乘法拟合最小二乘法拟合1 利用利用polyfit函数进行拟合函数进行拟合 polyfit函数的调用格式为：函数的调用格式为： P,S=polyfit(X,Y,m) 函数根据采样点函数根据采样点X和采样点函数值和

22、采样点函数值Y，产生一个，产生一个m次多项式次多项式P及其在采样点的误差向量及其在采样点的误差向量S。其中。其中X,Y是两是两个等长的向量，个等长的向量，P是一个长度为是一个长度为m+1的向量，的向量，P的元素的元素为多项式系数。为多项式系数。 例：例：x=0.5 1.0 1.5 2.0 2.5 3.0; y=1.75 2.45 3.81 4.80 8.00 8.60; a=polyfit(x,y,2) x1=0.5:0.05:3.0; y1=a(1)*x1.2 +a(2)*x1 + a(3); plot(x,y,*,x1,y1,-r) 2 利用矩阵除法解决复杂型函数的拟合利用矩阵除法解决复杂

23、型函数的拟合 例：用最小二乘法求一个形如例：用最小二乘法求一个形如y=a+bex的经的经 验公式验公式 ，其中，其中a、b为待定参数为待定参数 x=0:0.4:2.0; y=1.2 1.3 1.45 1.6 2.0 2.5; x1=exp(x) x1=ones(6,1),x1; %由由y=x1ab 得得 ab=x1y %其中其中 ab= x2=0:0.01:2.0; y2=ab(1)+ab(2)*exp(x2); plot(x,y,*,x2,y2,r) ab7.2 积分与微分积分与微分数值积分基本原理：数值积分基本原理： 求解定积分的数值方法多种多样，如简单的求解定积分的数值方法多种多样，如简

24、单的梯形法、辛普生梯形法、辛普生(Simpson) 法、牛顿柯特斯法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它法等都是经常采用的方法。它们的基本思想都是将整个积分区间们的基本思想都是将整个积分区间a,b分成分成n个个子区间子区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。这样求定积分问题就分解为求和问题。7.2.1 Newton-Cotes系列数值求积公式系列数值求积公式 若将积分区间若将积分区间a,b划分为划分为n等份等份,步长步长 h=(b-a)/n,选取等距节点选取等距节点 xk =a+kh,构造下构造下面的

25、插值型求公式：面的插值型求公式：In =(b-a)其中其中xk =a+kh, Ck(n) 由由Cotes系数表系数表7.6给出给出.当当n=0时时,即为矩形公式即为矩形公式(稍做变化稍做变化)当当n=1时时,即为梯形公式即为梯形公式当当n=2时时,即为即为Simpson公式公式当当n=3时时,即为即为Cotes公式公式sxfCknknk)(0)(1 矩形求积公式矩形求积公式 y=cumsum(x) 对对n维向量维向量x，返回返回n维向量维向量y，y(i)是是x的前的前i个分量之和。个分量之和。 例例:x=0 1 2;3 4 5;cumsum(x,1) %列求和列求和cumsum(x,2) %行

26、求和行求和2 梯形求积公式梯形求积公式 y=trapz(x) 对对n维向量维向量x，返回标量数值返回标量数值y，它的值是它的值是x的的各分量相隔为各分量相隔为1单位时各小梯形面积之和。单位时各小梯形面积之和。例：求积分例：求积分 e-0.5t sin(t+ /6)dt 解：建立被积函数解：建立被积函数 function y=fun1(t) y=exp(-0.5*t).*sin(t+pi/6); 03 d=pi/1000; syms t t=0:d:3*pi; y=exp(-0.5*t).*sin(t+pi/6); nt=length(t); int(y,0,3*pi) y=fun1(t); v

27、pa(ans,4) y=exp(-0.5*t).*sin(t+pi/6); sc=cumsum(y)*d; scf=sc(nt) z=trapz(y)*d3 自适应法自适应法Simpson法求积法求积 Q=quad(F,a,b) 求函数求函数F(x) 从从a到到b的相对误差为的相对误差为1e-3的积的积分近似值分近似值 I,n=quad(fname,a,b,tol,trace) 其中其中fname是被积函数名。是被积函数名。a和和b分别是定积分的下限和上限。分别是定积分的下限和上限。tol用来控制积分精度，缺用来控制积分精度，缺省时取省时取tol=0.001。trace控制是否展现积分过程，若

28、取非控制是否展现积分过程，若取非0则展则展现积分过程，取现积分过程，取0则不展现，缺省时取则不展现，缺省时取trace=0。返回参数。返回参数I即即定积分值，定积分值，n为被积函数的调用次数。为被积函数的调用次数。 例：例： dx function f=fun2(x)f=x./(4+x.2);quad(fun2,0,1) quad(x./(4+x.2),0,1)4 自适应法的自适应法的Cotes求积公式求积公式 quad8 自适应法的自适应法的Cote求积公式求积公式(高阶方法高阶方法)Q=quad8(F,a,b) 使用使用Newton-Cotes法则求自适应法则求自适应递归法求函数递归法求函

29、数F(X)从从a到到b的相对误差为的相对误差为1e-3的积分近的积分近似值似值,其中其中F为函数名为函数名.I,n=quad8(fname,a,b,tol,trace) 其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。 该函数可该函数可以更精确地求出定积分的值，且一般情况下函数调以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于用的步数明显小于quad函数，从而保证能以更高的函数，从而保证能以更高的效率求出所需的定积分值。效率求出所需的定积分值。10 xx2+4例：求积分例：求积分编制函数文件编制函数文件: fun.m func

30、tion f=fun(x) f=exp(-x/2);quad8(fun,1,3,1e-10)quad8(exp(-x/2),1,3,1e-10)3/ 21xed x例：例： 分别用分别用quad函数和函数和quad8函数求定积分的近似值，并在相同函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。的积分精度下，比较函数的调用次数。调用函数调用函数quad求定积分：求定积分：format long; fx=inline(exp(-x); I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766 n = 65 调用函数调用函数quad8求定积分：求定积分

31、：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)I = 0.28579444254754n = 33inline函数用于定义函数。函数用于定义函数。 inline函数的一般形式为：函数的一般形式为：FunctionName=inline(任何有效的任何有效的matlab表达式表达式, p1,p2 ,.) ,其中其中p1,p2 ,是出现在表达式中是出现在表达式中的所有变量名的所有变量名.求解求解F(x)=x2*cos(a*x)-b ,a,b是标量；是标量；x是向量是向量 输入输入: fx=inline(x .2*cos(a*x)-b

32、 , x,a,b); fx(pi/3 pi/3.5,4,1)输出为：输出为：g=-1.5483 -1.72597.2.2 二重定积分的数值求解二重定积分的数值求解 使用使用MATLAB提供的提供的dblquad函数就可以直接求函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：出上述二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace) 该函数求该函数求f(x,y)在在a,bc,d区域上的二重定积分。区域上的二重定积分。参数参数tol，trace的用法与函数的用法与函数quad完全相同。完全相同。例：例： 计算二重定积分计算二重定积分(1

33、) 建立一个函数文件建立一个函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于统计被积函数的调用次数用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2) 调用调用dblquad函数求解。函数求解。global ki;ki=0;I=dblquad(fxy,-2,2,-1,1) kiI = 1.57449318974494ki = 1038I=dblquad(exp(-x.2/2).*sin(x.2+y),-2,2,-1,1)7.2.6 微分和差分微分和差分1 数值微分和差分数值微分和差分在在MATLAB中，没

34、有直接提供求数值导数的函数，只中，没有直接提供求数值导数的函数，只有计算向前差分的函数有计算向前差分的函数diff，其调用格式为：，其调用格式为：diff(x) 计算向量计算向量X的向前差分的向前差分,其值为其值为 x(2)-x(1) x(3)-x(2) x(n)-x(n-1) diff(X) 对矩阵对矩阵X，计算矩阵计算矩阵X的向前差分，的向前差分，其值为矩阵其值为矩阵列的差分列的差分 X(2：n,:)-X(1:n-1,:) diff(X,n) 求求n阶的差分值，如果阶的差分值，如果n size(X,DIM),diff函数将先计算可能的连续差分值，直到函数将先计算可能的连续差分值，直到siz

35、e(X,DIM)=1. 然后然后diff 沿任意沿任意 n+1 维进行差分计算；维进行差分计算；即：即：计算计算X的的n阶向前差分。阶向前差分。例如：例如： diff(X,2)=diff(diff(X)。diff(X,n,DIM) 求求n阶的差分值，如果阶的差分值，如果n size(X,DIM),函函数返回空数组数返回空数组DX=diff(A,n,dim)：计算矩阵计算矩阵A的的n阶差分，阶差分，dim=1时时(缺缺省状态省状态)，按列计算差分；，按列计算差分；dim=2，按行计算差分。，按行计算差分。2 符号微分和差分符号微分和差分 diff(S) 对由对由findsym返回的自变量微分符号

36、表达式返回的自变量微分符号表达式S diff(S，V) 或或diff(S，sym（V）对自变量）对自变量v微分符微分符号表达式号表达式Sdiff(S，n) 对正整数对正整数n，函数微分函数微分S表达式表达式n次次diff(S,n,V) 和和diff(S,V , n）两种形式都可以被识别两种形式都可以被识别 例例 x=sym(x); diff(sin(x2)3 梯度函数梯度函数 gradient 近似梯度函数近似梯度函数fx,fy=gradient(F) 返回矩阵返回矩阵F的数值梯度，的数值梯度，fx相当于相当于 dF/dx, 为为x 方向的差分值。方向的差分值。fy相当于相当于dF/dy,为为

37、y方向的差分值。各个方向的点间隔方向的差分值。各个方向的点间隔设为设为1。当。当F为向量分值，为向量分值，DF=gradient(F) 为为一维梯度一维梯度7.3 求一元高次方程与线性方程组求一元高次方程与线性方程组7.3.1 直接解法直接解法1利用左除运算符的直接解法利用左除运算符的直接解法对于线性方程组对于线性方程组Ax=b，可以利用左除运算符，可以利用左除运算符“”求解：求解： x=Aban xn+an-1 xn-1+a0=0 p=an , an-1 , ,a0 x0=roots(p)n元线性方程组可表达为元线性方程组可表达为 AXb A为为nxn阶系数矩阵阶系数矩阵 b为为nx1阶常数

38、矩阵阶常数矩阵 X为为nx1阶变量矩阵阶变量矩阵 XAb 例例 求解下列方程组求解下列方程组2557. 03927. 02786. 04002. 01784. 04240. 00643. 03781. 01920. 03645. 01550. 01129. 04015. 03872. 02246. 04043. 02943. 03678. 01234. 04096. 04321432143214321xxxxxxxxxxxxxxxx解：在解：在MATLAB命令窗中输入：命令窗中输入： a=0.4096 0.1234 0.3678 0.2943 0.2246 0.3872 0.4015 0.11

39、29 0.3645 0.1920 0.3781 0.0643 0.1784 0.4002 0.2786 0.3927; b=0.4043 0.1550 0.4240 -0.2557; x=ab2利用矩阵的分解求解线性方程组利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有的矩阵分解有LU分解、分解、QR分解、分解、Cholesky分解，以及分解，以及Schur分解、分解、Hessenberg分解、分解、奇异分解等。奇异分解等。(1) LU分解分解矩阵的矩阵的

40、LU分解就是将一个矩阵表示为一个交换下三分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵中已经证明，只要方阵A是非奇异的，是非奇异的，LU分解总分解总是可以进行的。是可以进行的。MATLAB提供的提供的lu函数用于对矩阵进行函数用于对矩阵进行LU分解，其分解，其调用格式为：调用格式为：L,U=lu(X)：产生一个上三角阵：产生一个上三角阵U和一个变换形式和一个变换形式的下三角阵的下三角阵L(行交换行交换)，使之满足，使之满足X=LU。注意，。注意，这里的矩阵这里的矩阵X必须是方阵。必须是方阵。L,U,P

41、=lu(X)：产生一个上三角阵：产生一个上三角阵U和一个下三角和一个下三角阵阵L以及一个置换矩阵以及一个置换矩阵P，使之满足，使之满足PX=LU。当然。当然矩阵矩阵X同样必须是方阵。同样必须是方阵。实现实现LU分解后，线性方程组分解后，线性方程组Ax=b的解的解x=U(Lb)或或x=U(LP*b)，这样可以大大提高运算速度。，这样可以大大提高运算速度。例：例： 用用LU分解求解线性方程组。分解求解线性方程组。命令如下：命令如下：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分解的第分解的第

42、2种格式，命令如下：种格式，命令如下：L,U ,P=lu(A);x=U(LP*b)(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、一个、一个上三角矩阵

43、上三角矩阵R以及一个置换矩阵以及一个置换矩阵E，使之满，使之满足足XE=QR。实现实现QR分解后，线性方程组分解后，线性方程组Ax=b的解的解x=R(Qb)或或x=E(R(Qb)。例：例： 用用QR分解求解线性方程组。分解求解线性方程组。命令如下：命令如下：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(Qb)(3) Cholesky分解分解如果矩阵如果矩阵X是对称正定的，则是对称正定的，则Chole

44、sky分解将矩阵分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为上三角矩阵为R，则下三角矩阵为其转置，即，则下三角矩阵为其转置，即X=RR。MATLAB函数函数chol(X)用于对矩阵用于对矩阵X进行进行Cholesky分解，其调用格式为：分解，其调用格式为：R=chol(X)：产生一个上三角阵：产生一个上三角阵R，使，使RR=X。若。若X为非对称正定，则输出一个出错信息。为非对称正定，则输出一个出错信息。R,p=chol(X)：这个命令格式将不输出出错信息。：这个命令格式将不输出出错信息。当当X为对称正定的，则为对称正定的，则p=0

45、，R与上述格式得到的与上述格式得到的结果相同；否则结果相同；否则p为一个正整数。如果为一个正整数。如果X为满秩矩为满秩矩阵，则阵，则R为一个阶数为为一个阶数为q=p-1的上三角阵，且满足的上三角阵，且满足RR=X(1:q,1:q)。实现实现Cholesky分解后，线性方程组分解后，线性方程组Ax=b变成变成RRx=b，所以，所以x=R(Rb)。例：例： 用用Cholesky分解求解线性方程组。分解求解线性方程组。命令如下：命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;R=chol(A)A=4,-1,1; -1, 4.25, 2.

46、75;1, 2.75,3.5;b=13,6,0;R,p= chol(A)x=R(Rb)命令执行时，出现错误信息，说明命令执行时，出现错误信息，说明A为非正为非正定矩阵。定矩阵。 7.3.2 迭代解法的几种形式迭代解法的几种形式1 Jacobi迭代法迭代法对于线性方程组对于线性方程组Ax=b，如果，如果A为非奇异方阵，即为非奇异方阵，即aii0(i=1,2,n)，则可将，则可将A分解为分解为A=D-L-U，其，其中中D为对角阵，其元素为为对角阵，其元素为A的对角元素，的对角元素，L与与U为为A的下三角阵和上三角阵，于是的下三角阵和上三角阵，于是Ax=b化为：化为：x=D-1(L+U)x+D-1b

47、与之对应的迭代公式为：与之对应的迭代公式为：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例：例： 用用Jacobi迭代法求解下列线性方程组。设迭代迭代法求解下列线性方程组。设迭代初值为初值

48、为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)610272109103232121xxxxxxx2 Gauss-Seidel(G-S)迭代法迭代法将将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该式即为该式即为Gauss-Seid

49、el迭代公式。和迭代公式。和Jacobi迭迭代相比，代相比，Gauss-Seidel迭代用新分量代替旧迭代用新分量代替旧分量，精度会高些。分量，精度会高些。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例例: 用用Gauss-Seidel迭代法求解下列线性方程迭代法求解下列线性方程组。设迭代初值为组。设迭代初值为0，迭代精度为，迭代精度为10

50、-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)例例: 分别用分别用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)7.4 非线性方程数值求解非线性方程数值求解7.4.1 单变量非

51、线性方程求解单变量非线性方程求解 在在MATLAB中提供了一个中提供了一个fzero函数，可以函数，可以用来求单变量非线性方程的根。该函数的用来求单变量非线性方程的根。该函数的调用格式为：调用格式为： z=fzero(fname,x0,tol,trace)其中其中fname是待求根的函数文件名，是待求根的函数文件名，x0为搜索为搜索的起点。一个函数可能有多个根，但的起点。一个函数可能有多个根，但fzero函数只给出离函数只给出离x0最近的那个根。最近的那个根。tol控制结控制结果的相对精度，缺省时取果的相对精度，缺省时取tol=eps，trace 指定迭代信息是否在运算中显示，为指定迭代信息是

52、否在运算中显示，为1时显时显示，为示，为0时不显示，缺省时取时不显示，缺省时取trace=0。例例: 求求f(x)=x-10 x+2=0在在x0=0.5附近的根。附近的根。 步骤如下：步骤如下：(1) 建立函数文件建立函数文件funx.m。 function fx=funx(x) fx=x-10.x+2; (2) 调用调用fzero函数求根。函数求根。z=fzero(funx,0.5) z=fzero（x-10.x+2,0.5） z = 0.37587.4.2 非线性方程组的求解非线性方程组的求解 对于非线性方程组对于非线性方程组F(X)=0，用，用fsolve函数求其数值函数求其数值解。解。

53、fsolve函数的调用格式为：函数的调用格式为： X=fsolve(fun,X0,option)其中其中X为返回的解，为返回的解，fun是用于定义需求解的非线性是用于定义需求解的非线性方程组的函数文件名，方程组的函数文件名，X0是求根过程的初值，是求根过程的初值，option为最优化工具箱的选项设定。最优化工具为最优化工具箱的选项设定。最优化工具箱提供了箱提供了20多个选项，用户可以使用多个选项，用户可以使用optimset命命令将它们显示出来。如果想改变其中某个选项，令将它们显示出来。如果想改变其中某个选项，则可以调用则可以调用optimset()函数来完成。例如，函数来完成。例如，Disp

54、lay选项决定函数调用时中间结果的显示方式，其中选项决定函数调用时中间结果的显示方式，其中off为不显示，为不显示，iter表示每步都显示，表示每步都显示，final只显只显示最终结果。示最终结果。optimset(Display,off)将设定将设定Display选项为选项为off。例例: 求下列非线性方程组在求下列非线性方程组在(0.5,0.5) 附近的数值解。附近的数值解。 (1) 建立函数文件建立函数文件myfun.m。function q=myfun(p)x=p(1);y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y); (2) 在给定的初值在给