第六章-Matlab数值计算2

1、第六章 主要内容,数据处理与多项式计算 数值微积分 线性方程组求解 非线性方程与最优化问题求解 常微分方程的数值求解 稀疏矩阵,6.4 线性方程组求解,6.4.1 直接解法 利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运算符“”求解： x=Ab 例 用直接解法求解下列线性方程组。,A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; x=Ab %X=inv(A)*b,矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessen

2、berg分解、奇异分解等。 通过这些分解方法求解线性方程组的优点是速度快、可以节省存储空间。,2. 利用矩阵分解,矩阵的LU分解就是将一个矩阵A表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。,1. LU分解,L,U=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足A=LU。注意，这里的矩阵A必须是方阵。,例 对矩阵A用LU分解,命令如下： A=1 -1 1;5 -4 3;2 1 1 L,U=lu(A),结果如下： L = 0.2000 -0.0769 1.0000 1.0000 0 0 0.40

3、00 1.0000 0 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846,可以进一步去验证L*U,例 对矩阵A用LU分解,命令如下： A=1 -1 1;5 -4 3;2 1 1 L,U,P=lu(A),结果如下： L = 1.0000 0 0 0.4000 1.0000 0 0.2000 -0.0769 1.0000 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 P = 0 1 0 0 0 1 1 0 0,L,U,P=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换

4、矩阵P，使之满足PA=LU。当然矩阵A同样必须是方阵。,MATLAB提供的lu函数用于对矩阵进行LU分解： L,U=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足A=LU。注意，这里的矩阵A必须是方阵。 L,U,P=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。当然矩阵A同样必须是方阵。 实现LU分解后，线性方程组Ax=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。,例 用LU分解求解线性方程组。,A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; L

5、,U=lu(A); x=U(Lb) L,U,P=lu(A); x=U(LP*b),对矩阵A进行QR分解，就是把A分解为一个正交矩阵Q和 一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。,2. QR分解,Q,R=qr(A)：产生一个正交阵Q和一个上三角阵R，使之满足A=QR。注意，这里的矩阵A必须是方阵。,Q,R,E=qr(A)：产生一个正交矩阵Q、一个上三角阵R以及一个置换矩阵E，使之满足AE=QR。当然矩阵A同样必须是方阵。,例 对矩阵A用QR分解,命令如下： A=1 -1 1;5 -4 3;2 1 1 Q,R=qr(A),结果如下： Q = -0.1826 -0.1501 -0.971

6、7 -0.9129 -0.3412 0.2242 -0.3651 0.9279 -0.0747 R = -5.4772 3.4689 -3.2863 0 2.4427 -0.2456 0 0 -0.3737,例 对矩阵A用QR分解,命令如下： A=1 -1 1;5 -4 3;2 1 1 Q,R,E=qr(A),结果如下： Q = -0.1826 -0.1501 -0.9717 -0.9129 -0.3412 0.2242 -0.3651 0.9279 -0.0747 R = -5.4772 3.4689 -3.2863 0 2.4427 -0.2456 0 0 -0.3737 E = 1 0

7、0 0 1 0 0 0 1,QR分解 对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。 Q,R=qr(A)：产生一个一个正交矩阵Q和一个上三角矩阵R，使之满足A=QR。 Q,R,E=qr(A)：产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足AE=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) Q,R,E=qr

8、(A); x=E*(R(Qb),如果矩阵A是对称正定的，则Cholesky分解将矩阵A分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即A=RR。 R=chol(A)：产生一个上三角阵R，使RR=A。若A为非对称正定，则输出一个出错信息。 R,p=chol(A)：这个命令格式将不输出出错信息。当A为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。 实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x=R(Rb)。,3.Cholesky分解,例 设,A=2,1,1;1,2,-1;1,-1,3; R=chol(A) R*R R,

9、p=chol(A),例 用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),? Error using = chol Matrix must be positive definite. Error in = test at 3 R=chol(A),迭代解法非常适合求解大型系数矩阵的方程组。在数值 分析中，迭代解法主要包括 Jacobi迭代法、Gauss- Serdel迭 代法、超松弛迭代法和两步迭代法。,6.4.2 迭代解法,1. 迭代解法- Jacobi迭代法,将一组x代入右端，

10、可以立即得到另一组x。如果两组x相等，那么它就是方程组的解，不等时可以继续迭代。,对于线性方程组Ax=b，如果A为非奇异方阵，即aii0，(i=1,2,n)，则可将A分解为A=D-L-U； 其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：Dx=(L+U)x+b x=D-1(L+U)x+D-1b 与之对应的Jacobi迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b 如果序列x(k+1)收敛于x，则x必是方程Ax=b的解。,Jacobi迭代法的MATLAB函数文件Jacobi.m如下： function y,n=jacobi(A,b,x0,ep

11、s) if nargin=3 eps=1.0e-6; elseif nargin3 error return end D=diag(diag(A); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵,B=D(L+U); f=Db; y=B*x0+f; n=1; %迭代次数 while norm(y-x0)=eps x0=y; y=B*x0+f; n=n+1; end,例 分别用Jacobi迭代求解下列线性方程组，看是否收敛。,A=10,-1,0;-1,10,-2;0,-2,10; b=9,7,6; x,n=jacobi(A,b,0,

12、0,0,1.0e-6),2. Gauss Serdel迭代法,和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。,function y,n=gauserdel(A,b,x0,eps) if nargin=3 eps=1.0e-6; elseif nargin3 error return end D=diag(diag(A); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵,G=(D-L)U; f=(D-L)b; y=G*x0+f; n=1; %迭代次数 while norm(y-x0)=eps x

13、0=y; y=G*x0+f; n=n+1; end,例 用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。 在命令中调用函数文件gauserdel.m，命令如下： A=10,-1,0;-1,10,-2;0,-2,10; b=9,7,6; x,n=gauserdel(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=gauserdel(a,b,0;0;0),用Jacob

14、i迭代收敛，而Gauss-Serdel迭代不收敛。因此，使用迭代法时，要考虑算法收敛性。,线性方程组Ax=b的求解 当系数矩阵A满秩，方程有唯一解x=A-1b； 当b=0时为齐次方程组，总有零解，x=0为平凡解；当A的秩小于n，齐次方程有无穷多个非平凡解，用matlab中的函数null(A,r)可求基础解系； 当b0时 当rank(A)=rank(A,b)=n时，方程组有唯一解：x=Ab； 当rank(A)=rank(A,b)n时，方程组有无穷解，其通解=特解+齐次方程组的通解； 当rank(A)rank(A,b)时，方程组无解。,6.4.3 求线性方程组的通解,求线性方程组的函数文件 fun

15、ction x,y=line_solution(A,b) m,n=size(A); y=; if norm(b)0 if rank(A)=rank(A,b) if rank(A)=n disp(原方程组有唯一解x); x=Ab; else disp(原方程组有无穷个解，特解为x，其齐次方程组的基础解系为y); x=Ab; y=null(A,r); end else disp(方程组无解); x=; end else,disp(原方程组有零解x); x=zeros(n,1); if rank(A)n disp(方程组有无穷个解，基础解系为y); y=null(A,r); end end,A=1,

16、-2,3,-1;3,-1,5,-3;2,1,2,-2; b=1;2;3; x,y=line_solution(A,b) 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,y format short,6.5 非线性方程与最优化问题求解,6.5.1 单变量非线性方程求解 在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。该函数的调用格式为： z=fzero(fname,x0,tol,trace) 其中fname是待求根的函数文件名， x0为搜索的起点。一个函数可能有

17、多个根，但fzero函数只给出离x0最近的那个根； tol控制结果的相对精度，缺省时取tol=eps； trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。,例：分别用两个初始点求非线性方程的解,先建立函数文件fz function f=fz(x) f=x-1/x+5; 调用fz函数求根 fzero(fz,-5) fzero(fz,1) fzero(fz,3),ans = -5.1926 ans = 0.1926 ans = -5.1926,函数f 的图形 x=-10:0.1:5; y=x-1./x+5; plot(x,y),2. 非线性方程组的求解 对于

18、非线性方程组F(X)=0，用fsolve函数求其数值解。 X=fsolve(fun,X0,option) 其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值； option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。 例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。,例 求下列

19、非线性方程组在(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) 在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根。 x=fsolve(myfun,0.5,0.5,optimset(Display,off) x = 0.6354 0.3734,将求得的解代回原方程，可以检验结果是否正确，命令如下： q=myfun(x) q = 1.0e-009 * 0

20、.2375 0.2957 可见得到了较高精度的结果。,6.5.2 无约束最优化问题求解 无约束最优化问题的一般描述为minf(x)，寻求一组x，使得目标函数f(x)为最小，故又称最小化问题。 x,fval=fminbnd(filename,x1,x2,option):求一元函数在(x1,x2)区间中的最小值点x和最小值fval。 x,fval=fminsearch(filename,x0,option):基于单纯形算法求多元函数的最小值点x和最小值fval。 x,fval=fminunc(filename,x0,option):基于拟牛顿法求多元函数最小值点x和最小值fval。 x1和x2分别

21、表示研究区间的左右界；x0是一个向量，表示极值点的初值。,例 求函数在区间(-10,1)和(1,10)上的最小值点。 (1) 建立函数文件fx.m。 function f=fx(x) f=x-1/x+5; (2) 在Matlab命令窗口中输入命令： x,fmin=fminbnd(fx,-10,-1) x,fmin= fminbnd(fx,1,10),结果 x = -9.9999 fmin = -4.8999 ans = 1.0001,例 求函数在(0.5,0.5,0.5)附近上的最小值点。 (1) 建立函数文件fxyz.m。 function f=fxyz(u) x=u(1);y=u(2);z

22、=u(3); f=x+y.2./x/4+z.2./y+2./z; (2) 在Matlab命令窗口中输入命令： U,fmin=fminsearch(fxyz,0.5,0.5,0.5),结果 U = 0.5000 1.0000 1.0000 fmin = 4.0000,6.5.3 有约束最优化问题求解,约束条件可以进一步细化为： 线性不等式约束： 线性等式约束： 非线性不等式约束： 非线性等式约束： X的上下界：,x,fval=fmincon(filename,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,option),例 求解有约束的最优化问题。 (1) 建立函数文件fop.m。

23、 function f=fop(x) f=0.4*x(2)+x(1)2+x(2)2-x(1)*x(2)+1/30*x(1)3; (2) 编写函数： x0=0.5;0.5; A=-1,-0.5;-0.5,-1; b=-0.4;-0.5; lb=0;0; option=optimset;option.LargeScale=off;option.Display=off; x,f=fmincon(fop,x0,A,b,lb,option),结果 x = 0.3394 0.3303 f = 0.2456,龙格库塔法 基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为： t,y=o

24、de23(fname,tspan,y0) t,y=ode45(fname,tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量； tspan形式为t0,tf,表示求解区间； y0是初始状态列向量； t和y分别给出时间向量和相应的状态向量。,6.6 常微分方程的数值求解,例 求微分方程数值解，并与精确解相比较。 (1) 建立函数文件funt.m。 function y=funt(t,y) y=(y2-t-2)/4/(t+1); (2) 编写函数： t0=0;tf=10; y0=2; t,y=ode23(funt,t0,tf,y0); y1=sqrt(t+

25、1)+1; plot(t,y,b.,t,y1,r-);,例 求二阶微分方程数值解，并绘制时间响应曲线和相平面图。,(1) 建立函数文件funt.m function xdot=funt(t,x) xdot=-2*x(2);x(1); (2) 编写函数： t0=0;tf=20;x0=1,0; t,x=ode45(funt,t0,tf,x0); subplot(1,2,1);plot(t,x(:,2); subplot(1,2,2);plot(x(:,2),x(:,1),6.7 稀疏矩阵,(1) 矩阵存储方式 MATLAB的矩阵有两种存储方式：完全存储方式和稀疏存储方式。 完全存储方式完全存储方式

26、是将矩阵的全部元素按列存储。以前讲到的矩阵的存储方式都是按这个方式存储的，此存储方式对稀疏矩阵也适用。,稀疏存储方式 稀疏存储方式仅存储矩阵所有的非零元素的值及其位置，即行号和列号。在MATLAB中，稀疏存储方式也是按列存储的。 注意，在讲稀疏矩阵时，有两个不同的概念，一是指矩阵的0元素较多，该矩阵是一个具有稀疏特征的矩阵，二是指采用稀疏方式存储的矩阵。,函数A=sparse(S)将矩阵S转化为稀疏存储方式的矩阵A。当矩阵S是稀疏存储方式时，则函数调用相当于A=S； sparse(m,n)：生成一个mn的所有元素都是0的稀疏矩阵；,(2) 将完全存储方式转化为稀疏存储方式,程序： X=2 0

27、0 0 0;0 0 0 0 0;0 0 0 5 0;0 1 0 0 -1;0 0 0 0 5 A=sparse(X) 结果： A = (1,1) 2 (4,2) 1 (3,4) 5 (4,5) -1 (5,5) 5,例：将矩阵转化为稀疏存储方式,sparse(u,v,S)：u,v,S是3个等长的向量。S是要建立的稀疏矩阵的非0元素，u(i)、v(i)分别是S(i)的行和列下标；该函数建立一个max(u)行、max(v)列并以S为稀疏元素的稀疏矩阵； u,v,S=find(A)：返回矩阵A中非0元素的下标和元素。这里产生的u,v,S可作为sparse(u,v,S)的参数。 full(A)：返回和

28、稀疏存储矩阵A对应的完全存储方式矩阵。,(2) 将完全存储方式转化为稀疏存储方式,程序： A=sparse(1,2,5,7,1,3,2,3,1,1,2,4) B=full(A),例：生成稀疏矩阵和完全存储方式矩阵,结果： A = (1,1) 1 (5,2) 2 (2,3) 1 (7,3) 4 B = 1 0 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 4,只把要建立的稀疏矩阵的非0元素及其所在行和列的位置表示出来后由MATLAB自己产生其稀疏存储，调用格式为： B=spconvert(A) 其中A为一个m3或m4的矩阵，其每行表示一个非0元素，m是非0元素的个数，A每个元素的意义是：(i,1) 第i个非0元素所在的行。(i,2) 第i个非0元素所在的列。(i,3) 第i个非0元素值的实部。(i,4) 第i个非0元素值的虚部，若矩阵的全部元素都是实数，则无须第四列。 该函数将A所描述的一个稀疏矩阵转化为一个稀疏存储矩阵。,(3) 产生稀疏存储矩阵,例：根据表示稀疏矩阵的矩阵A，产生一个稀疏存储方式矩阵B。,命令如下：A=2,2,1;3,1,-1;4,3,3;5,3,8;6,6,12;B=spconvert(A) 结果如下： B = (3,1) -1 (2,2) 1 (4,3) 3 (5,3)