1、第 7 章 MATLAB 解方程与函数极值7.1 线性方程组求解7.2 非线性方程数值求解7.3 常微分方程初值问题的数值解法7.4 函数极值7.1 线性方程组求解7.1.1 直接解法1利用左除运算符的直接解法对于线性方程组 Ax=b，可以利用左除运算符“”求解：x=Ab例 7-1 用直接解法求解下列线性方程组。命令如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;x=Ab2利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有 LU 分解、QR 分解、 Cholesky 分解
2、，以及 Schur 分解、Hessenberg 分解、奇异分解等。(1) LU 分解矩阵的 LU 分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵 A 是非奇异的，LU 分解总是可以进行的。MATLAB 提供的 lu 函数用于对矩阵进行 LU 分解，其调用格式为：L,U=lu(X)：产生一个上三角阵 U 和一个变换形式的下三角阵 L(行交换) ，使之满足X=LU。注意，这里的矩阵 X 必须是方阵。L,U,P=lu(X)：产生一个上三角阵 U 和一个下三角阵 L 以及一个置换矩阵 P，使之满足PX=LU。当然矩阵 X 同样必须是方阵。实现 LU
3、分解后，线性方程组 Ax=b 的解 x=U(Lb)或 x=U(LPb)，这样可以大大提高运算速度。例 7-2 用 LU 分解求解例 7-1 中的线性方程组。命令如下：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 分解对矩阵 X 进行 QR 分解，就是把 X 分解为一个正交矩阵 Q 和一个上三角矩阵 R 的乘积形式。QR 分解只能对方阵进行。MATLAB 的函数 qr 可用于对矩阵进行 QR 分解，其调用
4、格式为：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)。例 7-3 用 QR 分解求解例 7-1 中的线性方程组。命令如下：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) C
5、holesky 分解如果矩阵 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
6、(1:q,1:q)。实现 Cholesky 分解后，线性方程组 Ax=b 变成 RRx=b，所以 x=R(Rb)。例 7-4 用 Cholesky 分解求解例 7-1 中的线性方程组。命令如下：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 = cholMatrix must be positive definite命令执行时，出现错误信息，说明 A 为非正定矩阵。7.1.2 迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Ser
7、del 迭代法、超松弛迭代法和两步迭代法。1Jacobi 迭代法对于线性方程组 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与之对应的迭代公式为：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=jaco
8、bi(A,b,x0,eps)if nargin=3eps=1.0e-6;elseif nargin=epsx0=y;y=B*x0+f;n=n+1;end例 7-5 用 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)2Gauss-Serdel 迭代法在 Jacobi 迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式 Dx(k+1)=(L+U)x(k)+b 可以改进为 Dx(k+1)
9、=Lx(k+1)+Ux(k)+b，于是得到：x(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=3eps=1.0e-6;elseif nargin=epsx0=y;y=G*x0+f;n=n+1;end例 7-6 用 Gauss-Serdel 迭代法求解下列线性方程组。设迭代初值为
10、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)例 7-7 分别用 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.2 非线性方程数值求解7.2.1 单变量非线性方程求解在 MATLAB 中提供了一个 fzero 函数，可以用来求单变
11、量非线性方程的根。该函数的调用格式为：z=fzero(fname,x0,tol,trace)其中 fname 是待求根的函数文件名，x0 为搜索的起点。一个函数可能有多个根，但 fzero函数只给出离 x0 最近的那个根。tol 控制结果的相对精度，缺省时取 tol=eps，trace 指定迭代信息是否在运算中显示，为 1 时显示，为 0 时不显示，缺省时取 trace=0。例 7-8 求 f(x)=x-10x+2=0 在 x0=0.5 附近的根。步骤如下：(1) 建立函数文件 funx.m。function fx=funx(x)fx=x-10.x+2;(2) 调用 fzero 函数求根。z=
12、fzero(funx,0.5)z =0.37587.2.2 非线性方程组的求解对于非线性方程组 F(X)=0，用 fsolve 函数求其数值解。fsolve 函数的调用格式为：X=fsolve(fun,X0,option)其中 X 为返回的解，fun 是用于定义需求解的非线性方程组的函数文件名， X0 是求根过程的初值，option 为最优化工具箱的选项设定。最优化工具箱提供了 20 多个选项，用户可以使用 optimset 命令将它们显示出来。如果想改变其中某个选项，则可以调用 optimset()函数来完成。例如，Display 选项决定函数调用时中间结果的显示方式，其中off为不显示，i
13、ter表示每步都显示， final只显示最终结果。optimset( Display,off)将设定Display 选项为off 。例 7-9 求下列非线性方程组在 (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
14、 =0.63540.3734将求得的解代回原方程，可以检验结果是否正确，命令如下：q=myfun(x)q =1.0e-009 *0.2375 0.2957可见得到了较高精度的结果。7.3 常微分方程初值问题的数值解法7.3.1 龙格库塔法简介7.3.2 龙格库塔法的实现基于龙格库塔法，MATLAB 提供了求常微分方程数值解的函数，一般调用格式为：t,y=ode23(fname,tspan,y0)t,y=ode45(fname,tspan,y0)其中 fname 是定义 f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan 形式为t0,tf,表示求解区间。y0 是初始状态列向量。t
15、 和 y 分别给出时间向量和相应的状态向量。例 7-10 设有初值问题，试求其数值解，并与精确解相比较( 精确解为 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; %求精确解tyy1y 为数值解，y1 为精确值，显然两者近似。例 7-11 求解著名的 Van der Pol 方程。例 7-12 有 Lorenz 模型的状态方程，试绘制系统相平面图。7.4 函数极值M
16、ATLAB 提供了基于单纯形算法求解函数极值的函数 fmin 和 fmins， 它们分别用于单变量函数和多变量函数的最小值，其调用格式为：x=fmin(fname,x1,x2)x=fmins(fname,x0)这两个函数的调用格式相似。其中 fmin 函数用于求单变量函数的最小值点。fname 是被最小化的目标函数名，x1 和 x2 限定自变量的取值范围。fmins 函数用于求多变量函数的最小值点，x0 是求解的初始值向量。MATLAB 没有专门提供求函数最大值的函数，但只要注意到-f(x)在区间(a,b)上的最小值就是 f(x)在 (a,b)的最大值，所以 fmin(f,x1,x2)返回函数 f(x)在区间 (x1,x2)上的最大值。例 7-13 求 f(x)=x3-2x-5 在0,5内的最小值点。(1) 建立函数文件 mymin.m。function fx=mymin(x)fx=x.3-2*x-5;(2) 调用 fmin 函数求最小值点。x=fmin(mymin,0,5)x=0.8165
