数学建模线性方程组的数值解法

大学数学实验MathematicalExperiments实验5线性代数方程组的数值解法

许多实际问题归结为线性（代数）方程组

大型（稀疏）方程组需要有效的数值解法

需要注意方程组的病态性问题数值解法(迭代解法)的收敛性为什么要学习线性方程组的数值解法机械设备、土建结构的受力分析输电网络、管道系统的参数计算经济计划企业管理3.线性方程组数值解法的MATLAB实现

4.实际问题中方程组的数值解1.两类数值解法:直接方法；迭代方法实验5的主要内容2.

超定线性方程组的最小二乘解线性方程组的一般形式、两类解法直接法经过有限次算术运算求出精确解（实际上由于有舍入误差只能得到近似解）----高斯（Gauss）消元法及与它密切相关的矩阵LU分解迭代法从初始解出发，根据设计好的步骤用逐次求出的近似解逼近精确解

----雅可比（Jacobi）迭代法和高斯—塞德尔（Gauss—Seidel）迭代法直接法-高斯消元法消元过程回代过程条件直接法-列主元消元法高斯消元法条件解决办法选最大的一个(列主元)将列主元所在行与第k行交换后,再按上面的高斯消元法进行下去，称为列主元消元法。(绝对值)很小时，用它作除数会导致舍入误差的很大增加直接法-高斯消元法的矩阵表示相当于方程AX=b两边左乘单位下三角阵M1高斯消元法的第一次消元直接法-高斯消元法的矩阵表示第二次消元相当于再左乘单位下三角阵M2直接法-矩阵LU分解高斯消元法通过左乘M,使MA=UM单位下三角阵，U上三角阵记L=M-1，L为单位下三角阵若A可逆且顺序主子式不为零，则A可分解为一个单位下三角阵L和一个上三角阵U的积A=LU。这种分解是唯一的，称矩阵LU分解。定理若A可逆，则存在交换阵P

使PA=LUL为单位下三角阵，U为上三角阵。第i行与第k行交换直接法-矩阵LU分解乘以初等交换阵P~交换阵（单位阵经若干次行交换）定理直接法-对称正定矩阵的分解直接法-三对角矩阵的LU分解在三次样条插值和其它一些计算中，会遇到求解系数矩阵A具有三对角形式的线性方程组，这时A的LU分解（假定分解存在）可表为:L和U的计算公式为线性方程组Ax=f可通过等价的两个三角形方程组Ly=f和Ux=y求解如下:直接法-三对角矩阵的LU分解追赶法方程组的病态性（1,1）01.201.121=+xxx22x120201.122121=+=+xxxxx对b的扰动敏感(2,0)向量和矩阵的范数度量向量、矩阵大小的数量指标向量范数矩阵范数范数-=2)(||||max2AAATl表示最大特征根向量和矩阵范数的相容性条件1）设b有扰动b，分析x的变化

xA的条件数越大,(由b的扰动引起的)x的变化可能越大条件数与病态性x的(相对)变化不超过b的(相对)变化的Cond(A)倍,也大致上是A的(相对)变化的Cond(A)倍。条件数与病态性2）设A有扰动A，分析x的变化

xA的条件数越大,(由A的扰动引起的)x的变化越大条件数大的矩阵是病态矩阵迭代法-一个例子迭代法-雅可比迭代迭代格式迭代法–

高斯-塞德尔迭代Jacobi迭代公式Gauss-Seideil迭代公式改进

迭代法的收敛性Jacobi迭代Gauss-Seideil迭代一般迭代形式（相容性）迭代法的收敛性其中||B||是任何一种矩阵范数迭代法-超松弛(SOR)迭代Gauss-Seideil迭代公式改进超松弛迭代低松弛迭代Gauss-Seideil迭代～迭代法-超松弛SOR迭代收敛充要条件若A对称正定SOR迭代------解大型稀疏矩阵方程组

超定线性代数方程组的最小二乘解实验1§2.1汽车刹车距离建立了刹车距离d与车速v的关系：方程个数超过了未知数个数，称为超定方程组

数据拟合已知一组数据，即平面上n个点(xi,yi),i=1,…n,

寻求一个函数（曲线）y=f(x),

使f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。

数据拟合问题的提法(#)+++++++++xyy=f(x)(xi,yi)i使点(xi,yi)与曲线y=f(x)的距离i尽量小,i=1,…n曲线拟合与最小二乘准则

先选定一组函数

r1(x),r2(x),…rm(x),m<n,令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…am

为待定系数。线性最小二乘法的基本思路记

按照最小二乘准则,问题归结为：求

a1,a2,…am

使J(a1,a2,…am)

最小。线性最小二乘法的求解记当RTR

可逆时(4)有唯一解线性最小二乘法的求解问题１怎样选择

{r1(x),

…rm(x)}，以保证系数{a1,…am}有唯一解？{a1,…am}有唯一解RTR可逆Rank(RTR)=mRank(R)=mR列满秩{r1(x),…rm(x)}在xi点线性无关(i=1,…n)分析为什么要规定m<n?若m=n或m>n,会如何?若m>n，a有无穷多解若m<n，超定方程组,a无解；求最小二乘解若m=n，R可逆时a有唯一解强令f(xi)=a1r1(xi)+…+amrm(xi)=

yi(i=1,…n)(即曲线f(x)过全部数据点，此时J=0)得问题２分析线性最小二乘中的“线性”指的是什么？问题３f(x)对a线性，于是求解线性方程组基函数{r1(x),…rm(x)}的选取1.通过机理分析建立数学模型来确定

f(x)++++++++++++++++++++++++++++++f=a1+a2x2.将数据(xi,yi)i=1,…n

作图，通过直观判断确定f(x)f=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=a1exp(a2x)f=a1exp(a2x)f

对a

非线性，怎么办？

线性方程组数值解法的MATLAB实现若A为可逆方阵，输出原方程的解x若A为nm矩阵（n>m），且ATA可逆，输出原方程的最小二乘解x

线性方程组数值解法的MATLAB实现[x,y,p]=lu(A)

若A可逆，输出x为单位下三角阵L，y为上三角阵U，p为一交换阵P，使PA=LU.u=chol(A)

对正定对称矩阵A的Cholesky分解，输出u为上三角阵U，使A=UTU2.矩阵LU分解

[x,y]=lu(A)

若A可逆且顺序主子式不为零，输出x为单位下三角阵L，y为上三角阵U，使A=LU；若A可逆，x为一交换阵与单位下三角阵之积.例.解A=[1031;2-103;1310],b=[14-514]',x=A\b,[L1,U1]=lu(A);L1,U1,A1=L1*U1,[L2,U2,P]=lu(A);L2,U2,P,A2=L2*U2,A3=inv(P)*A2并对系数矩阵作LU分解shiyan51若第1个方程改为3x2+x3=14结果如何当n很大时Hilbert矩阵呈病态

线性方程组数值解法的MATLAB实现3.范数条件数特征值

n=norm(x)

输入x为向量或矩阵，输出为x的2-范数

c=cond(x)

输入x为矩阵,输出为x的2-条件数

r=rcond(x)

输入x为方阵,输出为x条件数倒数

e=eig(x)

输入x为矩阵，输出x的全部特征值H=hilb(5),h=rats(H),b=ones(5,1);x=H\b;b(5)=1.1;x1=H\b;[x,x1],n1=cond(H),n2=rcond(H),例：观察Hilbert矩阵的病态性例.Hx=b,其中H=hilb(5),b=[1,…1]Tshiyan52xx11.0e+003*0.00500.0680-0.1200-1.38000.63006.3000-1.1200-9.94000.63005.0400cond(H)=4.7661e+0051.提取（产生）对角阵v=diag(x)

输入向量x，输出v是以x为对角元素的对角阵；输入矩阵x，输出v是x的对角元素构成的向量；例：v=diag(diag(x))

输入矩阵x，输出v是x的对角元素构成的对角阵，可用于迭代法中从A中提取D。2.提取上（下）三角阵

其他相关的MATLAB函数y=triu(x)

输入矩阵x，输出v是x的上三角阵；v=tril(x)

输入矩阵x，输出v是x的下三角阵；v=triu(x,1)

同上，但对角元素为0，可从A中提取U;v=tril(x,-1)

同上，但对角元素为0，可从A中提取L。例.用迭代法解shiyan53MATLAB对稀疏矩阵的处理:进行大规模计算的优点a=sparse(r,c,v,m,n)

在第r行、第c列输入数值v，矩阵共m行n列，输出a为稀疏矩阵，只给出(r,c)及v

aa=full(a)

输入稀疏矩阵a，输出aa为满矩阵（包含零元素）a=sparse(2,2:3,8,2,4),aa=full(a),a=(2,2)8aa=0000(2,3)80880输出n=500;b=[1:n]';a1=sparse(1:n,1:n,4,n,n);a2=sparse(2:n,1:n-1,1,n,n);a=a1+a2+a2';tic;x=a\b;t1=tocaa=full(a);tic;xx=aa\b;t2=tocy=sum(x)yy=sum(xx)例.分别用稀疏矩阵和满矩阵求解Ax=b,比较计算时间设00t1,t2相差巨大，说明用稀疏矩阵计算的优点（y=yy用于简单地验证两种方法结果的一致）shiyan54用MATLAB作线性最小二乘拟合1.作多项式f(x)=a1xm+…+amx+am+1拟合,可利用已有程序:a=polyfit(x,y,m)输入:数据x,y(同长度数组);m(拟合多项式次数)输出:系数a=[a1,…am,

am+1](数组)。2.对超定方程组仍用可得最小二乘意义下的解多项式在x点的值:

y=polyval(a,x)实例1投入产出模型表2投入产出表假定每个部门的产出与各部门对它的投入成正比，得到投入系数。4）如果对于任意给定的、非负的外部需求，都能得到非负的总产出，模型就称为可行的。问为使模型可行，投入系数应满足什么条件？1）设有n个部门，已知投入系数，给定外部需求，建立求解各部门总产出的模型。2）设投入系数如表2所给，如果今年对农业、制造业和服务业的外部需求分别为50，150，100亿元，问这三个部门的总产出分别应为多少。3）如果三个部门的外部需求分别增加1个单位，它们的总产出应分别增加多少。实例1投入产出模型1）基本模型xi:第i个部门的产出，xij:第i个部门对第j个部门的投入，di:第i个部门的外部需求投入系数基本模型2）设农业、制造业和服务业的外部需求分别为50，150，100亿元，求三个部门的总产出。x=(139.2801,267.6056,208.1377)Tshiyan55若d=(1,0,0)T,

即农业外部需求增加1单位时，三部门总产出应分别增加1.3459，0.5634，0.4382单位。即C的第1列。C的第2，3列给出了什么？C=1.34590.25040.34430.56341.26760.49300.43820.43041.21673）若三部门的外部需求分别增加1个单位，求它们的总产出的增量。基本模型当需求增加d时，总产出增量4）如果对于任意给定的、非负的外部需求，都能得到非负的总产出，模型就称为可行的。问为使模型可行，投入系数应满足什么条件？模型可行模型可行若A是由实际数据得到，只要初始投入非负模型可行实验