计算机应用基础-3-线性与非线性方程组求解_精品文档.ppt

1、第三章 代数方程（组）的Matlab求解,3.1 线性方程（组）求解一方程组解个数的判定二高斯消元法求解三LU矩阵分解求解四Jacobi迭代求解五 调用内部函数求解3.2 非线性方程（组）求解,3.1 线性方程（组）求解,一、方程组解个数及求解,Ax=B,对于矩阵,3.1 线性方程（组）求解,解的判定矩阵,又叫增广矩阵,线性方程组解的判定定理:(1)当 m=n 且 rank(A)=rank(C)=n 时 方程组有唯一解 x=A-1B x=inv(A)*B(2)当rank(A)=rank(C)=rn 时,方程组有无穷多解:,Z=null(A)%求解A矩阵的化零矩阵Z=null(A,r)%求解A矩

2、阵的化零矩阵的规范形式,若rank(A)=rank(C),只能用Moore-penrose 广义逆求解方程的最小二乘解为:x=pinv(A)*B 只能使误差范数测度 取最小值。,3.1 线性方程（组）求解,【例3-1】求解,X=-1.8000 2.4000 1.8667-1.2667 3.8667-3.2667-2.1333 2.7333,3.1 线性方程（组）求解,检验：norm(A*x-B)=7.47e-015,1 2 3 4 5 1 4 3 2 1 X=4 2 1 3 2 4 3 3 4 1 3 2 2 4,A=1 2 3 4;4 3 2 1;1 3 2 4;4 1 3 2;B=5 1;

3、4 2;3 3;2 4;X=inv(A)*B,X1=-9/5,12/5 28/15,-19/15 58/15,-49/15-32/15,41/15,3.1 线性方程（组）求解,精确解：x1=inv(sym(A)*B,检验：norm(double(A*x1-B)=0,【例3-2】,A=1 2 3 4;2 2 1 1;2 4 6 8;4 4 2 2 B=1;3;2;6;C=A B;rank(A),rank(C),rank(A)=rank(C)=24,3.1 线性方程（组）求解,1 2 3 4 1 2 2 1 1 X=3 2 4 6 8 2 4 4 2 2 6,判定可解性：,Z=null(A,r),

4、Z=2.0000 3.0000-2.5000-3.5000 1.0000 0 0 1.0000,X0=0.9542;0.7328;-0.0763;-0.2977,X0=pinv(A)*B,3.1 线性方程（组）求解,解出规范化零空间：,得出一个特解：,3.1 线性方程（组）求解,解出全部解：,2 3 0.9542X=a1-2.5+a2-3.5+0.7328 1 0-0.0763 0 1-0.2977,验证得出的解：a1=randn(1);a2=rand(1);x=a1*Z(:,1)+a2*Z(:,2)+x0 norm(A*x-B),【例3-3】,3.1 线性方程（组）求解,1 2 3 4 1

5、2 2 1 1 X=2 2 4 6 8 3 4 4 2 2 4,B=1:4;C=A B;rank(A),rank(C),最小二乘解：x=pinv(A)*B,检验：norm(A*x-B)=7.47e-015,二 高斯消去法求解线性方程组 高斯消去法它不需要方程组解的初值，也不需要重复迭代计算，它通过消去和回代两个过程就可以直接求出方程组的解。和前面的迭代法一样，在计算过程中，也有可能发散而得不到方程组的解，但可以通过一些其它方法解决。设一个n元方程组，以矩阵形式表示为：,3.1 线性方程（组）求解,1 消元过程,（1）令，进行矩阵的代数运算，可消去第一列的元素（2 n),(2)同理，可消去第二列

6、的元素（3 n),（3）可得任意步消元公式，假设第k步，有：,3.1 线性方程（组）求解,（4）经过n-1步消元后，系数矩阵称为上三角矩阵,3.1 线性方程（组）求解,回归过程：（1）首先解最后一个方程（2）依次回代，可得到任意公式:,3.1 线性方程（组）求解,高斯消元法通用程序,function x=UpTriangleFun(A,b)%求上三角系数矩阵的线性方程组Ax=b的解%A为线性方程组的系数矩阵%b为线性方程组中的常数向量%x为线性方程组中的解n=size(A);N=n(1);for k=N:-1:1 if(kN)s=A(k,(k+1):N)*x(k+1):N,1);else s=

7、0;end x(k,1)=(b(k)-s)/A(k,k);end,3.1 线性方程（组）求解,function x=DownTriangleFun(A,b)%求下三角系数矩阵的线性方程组Ax=b的解%A为线性方程组的系数矩阵%b为线性方程组中的常数向量%x为线性方程组中的解n=size(A);N=n(1);for k=1:N if(k1)s=A(k,1:(k-1)*x(1:(k-1),1);else s=0;end x(k,1)=(b(k)-s)/A(k,k);end,3.1 线性方程（组）求解,高斯全主元消去法,对于求解精度较高，或方程数量较多，易采用全主元法,其增广矩阵为(A(1)|b(1

8、),寻找A(1)中绝对值最大为主元素,将此元素进行行列交换，换为第一个元素，再进行消元过程,3.1 线性方程（组）求解,对于任意的第k-1次，其最大的主元素为：,进过n-1次主元消去之后，对后得到上三角矩阵。,通用算法,3.1 线性方程（组）求解,function x,xA=AllGaussFun(A,b)%AllGaussFun.m为用高斯全主元消去法求解线性方程组的解%A为线性方程组的系数矩阵%b为线性方程组的常数向量%x为线性方程组的解%xA为消元后的系数矩阵N=size(A);n=N(1);indexl=0;indexr=0;order=1:n;%记录未知数顺序的向量for i=1:(

9、n-1)me=max(max(abs(A(i:n,i:n);%选取全主元 for k=i:n for r=i:n if(abs(A(k,r)=me)indexl=k;indexr=r;k=n;break;end end end,temp=A(i,1:n);A(i,1:n)=A(indexl,1:n);A(indexl,1:n)=temp;bb=b(indexl);b(indexl)=b(i);b(i)=bb;%交换主行 temp=A(1:n,i);A(1:n,i)=A(1:n,indexr);A(1:n,indexr)=temp;%交换主列 pos=order(i);order(i)=orde

10、r(indexr);order(indexr)=pos;%主列的交换会造成未知数顺序的变化 for j=(i+1):n if(A(i,i)=0)disp(对角元素为0！);return;end l=A(j,i);m=A(i,i);A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;b(j)=b(j)-l*b(i)/m;endend,x=UpTriangleFun(A,b);y=zeros(n,1);for i=1:n for j=1:n if(order(j)=i)y(i)=x(j);end endend x=y;xA=A;,例【3-4】求多元线性方程组的解,求解程序“li3_4.m

11、”,x=（2.40,-19.20,-1.00,4.00),3.1 线性方程（组）求解,三 LU矩阵分解求解,线性方程组变为：,LUx=b,令：,Y=Ux,则：,LY=b,3.1 线性方程（组）求解,解两个三角行方程组，可得解。,A=LU,其中,function x=LUfun(A,b)%实现LU分解求线性方程组Ax=b的解%A为线性方程组的系数矩阵%b为线性方程组中的常数向量%x为线性方程组中的解r=rank(A);%flag=isexist(A,b)%判断方程组解的情况if flag=0 disp(该方程组无解！);x=;return;else m,n=size(A);L,U,P=lu(A)

12、;%lu函数为MATLAB提供的函数实现矩阵LU分解 b=P*b;%解Ly=b y(1)=b(1);if m1 for k=2:m y(k)=b(k)-L(k,1:k-1)*y(1:k-1);end end y=y;,%解Ux=y得原方程组的一个特解 x0(r)=y(r)/U(r,r);if r1 for k=r-1:-1:1;x0(k)=(y(k)-U(k,k+1:r)*x0(k+1:r)/U(k,k);end end x0=x0;%如果方程组有唯一解 if flag=1 x=x0;return;else%如果方程组有无穷多解 format rat;z=null(A,r);%求出对应齐次方程

13、组的基础解系 mz,nz=size(z);x0(r+1:n)=0;for k=1:nz t=sym(char(107 48+k);K(k)=t;%取K=k1,k2,.;end x=x0;for k=1:nz x=x+K(k)*z(:,k);%将方程组的通解表示为特解加对应齐次通解形式 end endend,例 利用LU分解求解如下线性方程组的解,3.1 线性方程（组）求解,四 Jacobi迭代法求解,迭代法是对任意给定的初始向量，按照某种方法逐步生成近似解序列，是解序列的极限为方程组的解。,设线性方程组：Ax=b,若A非奇异，b不等于0，有唯一解x*，上式变换为：,x=Tx+C,取一初始向量x

14、(0)为上述方程的近似解,x（k+1）=Tx(k)+C,满足条件：,3.1 线性方程（组）求解,此迭代法收敛，x*为原方程组的解。,Jacobi迭代法 工程常用此方法解决稀疏矩阵的迭代问题。,系数矩阵为,3.1 线性方程（组）求解,若A非奇异，方程组Ax=b有唯一解x=A-1b.方程组通过变换为：,令,3.1 线性方程（组）求解,令,有：,3.1 线性方程（组）求解,算法：（1）设定一个初始向量x(0),将初始值带入矩阵方程，得到新向量x(1）。（2）将新向量x（1）带入矩阵方程，得到另一个向量x(2)。(3)对于任意k步迭代，迭代格式为：,如果序列x(k)收敛于x*，则说明x*为方程组的解,

15、Jacobi迭代程序：,3.1 线性方程（组）求解,function x,n=Jacobifun(A,b,x0,eps,var)%Jacobifun为编写的雅可比迭代函数%A为线性方程组的系数矩阵%b为线性方程组的常数向量%x0为迭代初始向量%eps为解的精度%var为迭代步数控制%x为线性方程组的解%n为求出所需精度的解实际的迭代步数if nargin=3 eps=1.0e-6;M=200;elseif nargin3 error returnelseif nargin=5 M=var1;end,3.1 线性方程（组）求解,D=diag(diag(A);L=-tril(A,-1);U=-tr

16、iu(A,1);B=D(L+U);f=Db;x=B*x0+f;n=1;%统计迭代次数while norm(x-x0)=eps x0=x;x=B*x0+f;n=n+1;if(n=M)disp(Warning:迭代次数太多，不收敛！);return;endend,3.1 线性方程（组）求解,例【3-5】用Jacobi迭代法计算方程组,设初始值为x0=0,0,0计算程序“li3_5.m”X=2.39,2.18,-0.25,3.1 线性方程（组）求解,五 调用Matlab内部函数,5.1 共轭梯度法，通过最优化求解，适用于系数矩阵为方阵m,调用格式如下：,x,flag,relres,iter=bicg(A,b,tol),方程组的解,方程组成功的标志,解的相对误差,迭代的次数,求解器,系数矩阵,误差限,常数向量,3.1 线性方程（组）求解,A=4,1,-1;2,-3,1;1,2,-5;b=12,-2,8;tol=1e-7;x,flag,relres,iter=bicg(A,b,tol),x=2.3929，2.1786-0.2500flag=0 relres=1.1384e-014 iter=3,