1、MATLAB的方程(组)解法,第六讲,王文健 ,MATLAB数据处理与应用2011-2012学年选修课,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,2,主要内容,线性方程(组)的解法 非线性方程(组)解法 MATLAB统计分析,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,3,简介,线性方程(组) 直接法在没有舍入的情况下，通过有限步四则运算求得方程组的准确解 迭代法先给出一个解的初始值，然后按一定的法则逐步求出解的各个更准确的近似值的方法 非线性方
2、程(组) 迭代法不动点迭代法、Newton迭代法、Broyden Broyden迭代法,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,4,线性方程(组)的解法,线性方程的解法 对于线性方程可以直接调用roots函数来求解 例如：求解x3+5x2-8x+6=0 P=1 5 -8 6; roots(P) 求解2x9+43x7+x6-8x5+14x4-5x3+x2-10x+12=0 a=2 0 43 1 -8 14 -5 1 -10 12; roots(a),Tribology Research Institute SOUTHW
3、EST JIAOTONG UNIVERSITY,5,线性方程(组)的解法,线性方程组的解法 直接法利用符号运算“/”或“”来求解 直接法一般基于高斯消去法、主元素消去法、平方根法和追赶法等，在MATLAB中这些算法已被编成现成的库函数或运算符，可以直接调用进行求解,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,6,线性方程(组)的解法,线性方程组的解法 建立系数矩阵和常数矩阵 A=2 1 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6 B=8 9 -5 0 X=AB,Tribology Research
4、Institute SOUTHWEST JIAOTONG UNIVERSITY,7,线性方程(组)的解法,线性方程组的解法 利用矩阵的LU、QR和Cholesky分解法求解这三种方法对求解大型方程组非常有用 优点是运算速度快、可以节省磁盘空间、节省内存 LU分解法Gauss消去分解，可以把任意矩阵分解为下三角矩阵的基本变换形式和上三角矩阵的乘积 A=LU L为下三角矩阵，U为上三角矩阵 A*X=b变成L*U*X=bX=U(Lb),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,8,线性方程(组)的解法,LU分解法 A=2 1
5、 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6 B=8 9 -5 0 L,U=lu(A) X=U(LB),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,9,线性方程(组)的解法,Cholesky分解法 如果A为对称正定矩阵，则Cholesky分解可将矩阵A分解为上三角矩阵和其转置的乘积，即A=R*R，其中R上三角矩阵 A*X=b变成R*R*X=b X=R(Rb) 举例： A=16 4 8;4 5 -4;8 -4 22 B=28 5 26 R=chol(A) X=R(RB),Tribology Resea
6、rch Institute SOUTHWEST JIAOTONG UNIVERSITY,10,线性方程(组)的解法,QR分解法 对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵，即A=QR A*X=b变成Q*R*X=b X=R(Qb) 举例： A=16 4 8;4 5 -4;8 -4 22 B=28 5 26 Q,R=qr(A) X=R(QB),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,11,线性方程(组)的解法,迭代解法 在MATLAB中迭代法非常适合求解大型系数矩阵的方程组 Jacobi迭代
7、法 Gauss-Serdel迭代法 超松弛迭代法 两步迭代法,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,12,线性方程(组)的解法,Jacobi迭代法 线性方程组Ax=b，如果系数矩阵A为非奇异矩阵，则A可写成A=D.L.U的分解形式，其中D为对角矩阵，元素为A的对角元素，L和U分别为严格的下三角矩阵和上三角矩阵，故迭代形式为： xi+1=D-1(L+U)xi+D-1b 对应的Jacobi迭代公式的矩阵为： x(k+1)=Bx(k)+f 式中：B=D-1(L+U)=I-D-1A f=D-1b,Tribology Res
8、earch Institute SOUTHWEST JIAOTONG UNIVERSITY,13,线性方程(组)的解法,Jacobi迭代法 举例：A=10,-1,0;-1 10 -2;0 -2 10; b=9;7;6; jacobi(A,b,0;0;1),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,14,线性方程(组)的解法,Gauss-Seidel迭代法 线性方程组Ax=b，Gauss-Seidel迭代形式为： X(k+1)=GX(k)+f 式中：G为Gauss-Seidel迭代矩阵，G=-(D+L)-1b，f= =(
9、D+L)-1b，D为对角矩阵，L和U为严格的下三角和上三角矩阵 编制Gauss-Seidel迭代法的M文件,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,15,线性方程(组)的解法,Gauss-Seidel迭代法 举例：A=10,-1,0;-1 10 -2;0 -2 10; b=9;7;6; GS(A,b,0;0;0),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,16,线性方程(组)的解法,SOR迭代法 线性方程组Ax=b，SOR迭代形式为： X(k
10、+1)=B0X(k)+f 式中：B0为SOR迭代矩阵，B0=(D-wL)-1(1-w)D+wU, f=w(D-wL)-1b D为对角矩阵，L和U为严格的下三角和上三角矩阵 编制SOR迭代法的M文件,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,17,线性方程(组)的解法,SOR迭代法 举例：A=4 3 -1;3 4 -1;1 -1 4; b=19;30;27; w=1.03; SOR(A,b,1;1;1,w),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSIT
11、Y,18,线性方程(组)的解法,双步迭代法 线性方程组Ax=b，如果A为对称正定矩阵，双步迭代法形式为： (D-L)x(k+1/2)=Ux(k)+b (D-U)x(k+1)=Lx(k+1/2)+b 式中：D为对角阵，L和U为严格的下三角和上三角矩阵 编制SOR迭代法的M文件,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,19,线性方程(组)的解法,双步迭代法 举例：A=10 -1 2 2;-1 11 -1 3;2 -1 10 3;2 3 -1 8; b=8;25;-11;17; TS_D(A,b,0;0;0;0),Trib
12、ology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,20,非线性方程的解法,对比法 假定函数在a,b区间上是单调函数，且f(a)f(b)0，则在区间上方程f(x)=0至少有一个根 思想：通过判断函数f(x)的符号，逐步将有限区间缩小，使在足够小的区间内方程有且仅有一个根，近似根就是最终区间的重点位置 举例：求解非线性方程xex-2=0 1.编制函数文件exam.m 2.命令窗口输入：DF(exam,0,10,0.05),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,
13、21,非线性方程的解法,迭代法 不定点迭代法 思想：非线性方程f(x)=0，如果在区间上连续，可将函数表示为x=g(x)，迭代格式为：xk+1=g(xk)，如果limxk=x*，则迭代过程将收敛于方程的根 举例：求解非线性方程ex-3x2+3=0 1.编制函数文件fun1.m 2.命令窗口输入：x,y,n=iterate(10),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,22,非线性方程的解法,迭代法 切线迭代法牛顿迭代法 思想：迭代格式xk+1=xk-f(xk)/fg(xk)-给出导数 举例：求解非线性方程ex-3x
14、2=0 1.编制函数文件ff.m 2.编制导数文件dff.m 3.命令窗口输入：x,n=NT(2),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,23,非线性方程的解法,迭代法 割线迭代法 思想：当不存在解析的导函数时，割线法可近似来求解导函数，迭代格式xk+1=xk-f(xk)(xk-xk-1)/(f(xk)- f(xk-1) 举例：求解非线性方程ex-3x2=0 1.编制函数文件GX.m 2.编制函数文件ff.m 3.命令窗口输入：y=GX(2,1),Tribology Research Institute SOUTH
15、WEST JIAOTONG UNIVERSITY,24,非线性方程的解法,求解非线性方程的MATALB函数 求解非线性方程f(x)=0的函数fsolve调用格式： x=fsolve(fun,x0)变量fun为方程表达式或函数文件名，x0为估计的方程根的初始值 x=fsolve(fun,x0,options)根据优化参数options进行最小化，默认值为0，如为1显示运算的中间步骤 x,fval=fsolve(fun,x0) 返回值为目标函数在求得的根x对应的函数值,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,25,非线性
16、方程的解法,求解非线性方程的MATALB函数 求解非线性方程f(x)=0的函数fsolve调用格式： x,fval,exitflag=fsolve()返回值为函数调用的退出标记exitflag，用于描述当前退出状态 x,fval,exitflag,output=fsolve()函数返回值为一个结构输出，用于描述优化的信息 x,fval,exitflag,output,jacobian=fsolve()函数返回值为在点x处函数的Jacobian数值,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,26,非线性方程的解法,求解非
17、线性方程的MATALB函数 举例：求解f(x)=xex-1=0的数值解 命令窗口输入： xx=0.5; options(1)=1;%优化过程显示运算的中间步骤 xx=fsolve(x*exp(x-1),xx,options)x,fval,exitflag,output,Jacobian=fsolve(x*exp(x-1),xx),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,27,非线性方程的解法,求解非线性方程的MATALB函数 fzero函数调用格式： x=fzero(fun,x0)fun为求解的方程，x0为估计的方程
18、的根或根可能存在的区间。 如果x0为一个数值，此函数就在x0附近求出方程的根，如果x0为一个区间，且函数的根包含在该区间内则输出正常的方程根，否则输出错误信息 x=fzero(fun,x0,options)通过options指定的优化参数进行求解 x,fval=fzero(fun,x0) fval该函数返回值为方程根处的函数值,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,28,非线性方程的解法,求解非线性方程的MATALB函数 fzero函数调用格式： x,fval,exitflag=fzero()exitflag返回函
19、数fzero退出的函数值 x,fval,exitflag,output=fzero()output返回优化的信息,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,29,非线性方程的解法,求解非线性方程的MATALB函数 举例：求解f(x)=3x2-exp(x)+1=0的根 命令窗口输入： x0=4; 赋初始值 fzero(3*x2-exp(x)+1,x0) x1=4 10; 区间未包含根 fzero(3*x2-exp(x)+1,x1) x1=3 5; 区间包含根 fzero(3*x2-exp(x)+1,x1),Tribolo
20、gy Research Institute SOUTHWEST JIAOTONG UNIVERSITY,30,非线性方程组的解法,非线性方程组的求解 求解更为复杂、繁琐，一般要采用迭代法进行求解 不动点迭代法 牛顿迭代法 Broyden迭代法,Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,31,非线性方程组的解法,不动点迭代法 思想：有n个未知数和n个方程的非线性方程组，F(x)=0，其迭代形式为：x=g(x)，不动点的迭代公式 x(k+1)=g(xk) 如果得到的序列x(k)满足limx(k)=x*，则x*为函数g(xk
21、)的不动点 举例：求解非线性方程组x1-0.6sinx1-0.3cosx2=0和x2-0.6cosx1+0.3sinx2=0 1.编制文件BDD.m 2.编制函数文件gg.m 3.命令窗口输入：y,n=BDD(0.5 0.5),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,32,非线性方程组的解法,牛顿迭代法 思想：有n个未知数和n个方程的非线性方程组，F(x)=0，其迭代形式为：x=g(x)，不动点的迭代公式 x(k+1)=x(k)-(F(x(k)-1F(x(k) 牛顿迭代法就是需要求解出线性方程组F(x(k)x(k)=
22、- F(x(k)的x(k)，如果x(k)足够小，则得到x(k+1)即为非线性方程组的解 举例：求解非线性方程组x12-10x1+x22+8=0和x1x22+x1-10x2+8=0 1.编制文件NewDD.m 2.编制函数文件fx1.m、dfx1.m 3.命令窗口输入：y,n=NewDD(0 0),Tribology Research Institute SOUTHWEST JIAOTONG UNIVERSITY,33,非线性方程组的解法,Broyden迭代法 思想：牛顿迭代法具有较好的收敛性，但每步都要计算Fg(x(k)的值，用较简单的矩阵A(k)年来代替牛顿迭代法中的Fg(x(k) ，就可以较为方便的进行迭代计算 举例：求解非线性方程组x1-0.7sinx1-0.2cosx2=0和x2-0.7cosx1+0.2sinx2=0 1.编制文件Broy.m 2.编制函数文件bf.m 3.命令窗口输入：x0=0.5 0.5; y,n=Broy(x0),谢 谢,
