1 大学数学实验 Experiments in Mathematics 实验5 线性代数方程组的数值解法 清华大学数学科学系 ? 许多实际问题归结为线性（代数）方程组 ? 大型的方程组需要有效的数值解法 ? 数值解法的稳定性和收敛性问题需要注意 为什么要学习线性方程组 的数值解法 机械设备、土建结构的受力分析 输电网络、管道系统的参数计算 经济计划 企业管理 3. 线性方程组数值解法的MATLAB实现 4. 实际问题中方程组的数值解。 1. 两类数值解法: 直接方法；迭代方法 实验 5的主要内容 2. 超定线性方程组的最小二乘解 线性方程组的一般形式、两类解法 直接法 经过有限次算术运算求出精确解（实际上 由于有舍入误差只能得到近似解）---- 高斯 （Gauss）消元法及与它密切相关的矩阵LU分解 迭代法 从初始解出发，根据设计好的步骤用逐次 求出的近似解逼近精确解 ---- 雅可比（Jacobi） 迭代法和高斯 —塞德尔（Gauss —Seidel）迭代法 nnnnnn nn nn bxaxaxa bxaxaxa bxaxaxa =+++ =+++ =+++ L LL L L 2211 22222121 11212111 或 AX=b 直接法---高斯消元法 )2()2( 2 )2( 2 )2( 2 )2( 22 )2( 22 )1( 1 )1( 12 )1( 121 )1( 11 nnnnn nn nn bxaxa bxaxa bxaxaxa =++ =++ =+++ L LL L L )()( )1( 1 )1( ,11 )1( 1,1 )2( 2 )2( 22 )2( 22 )1( 1 )1( 12 )1( 121 )1( 11 n nn n nn n nn n nnn n nn nn nn bxa bxaxa bxaxa bxaxaxa = =+ =++ =+++ ? ? ? ?? ? ?? LL LL LL 消 元 过 程 回 代 过 程 ),,2,1( 0 )( nk a k kk L= ≠ 条件 直接法-列主元素消元法 ),,2,1(0 )( nka k kk L=≠ 高斯消元法 条件 解决 办法 ),( )( nkia k ik L=选 最大的一个（列主元） 将列主元所在行与第 k行交换后 , 再按上面的高斯消元 法进行下去 ， 称为 列主元素消元法 。 (绝对值 )很小时， 用它作除数会导致舍入误 差的很大增加 )(k kk a )()()( )()()( k nn k nnk k nk k k n k kn k k kk bxaxa bxaxa =++ =++ LL LL LL 2 直接法 - 高斯消元法的矩阵表示 相当于方程 AX=b 两边 左乘单位下 三角阵 M 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = 1/ 1/ 1 )1( 11 )1( 1 )1( 11 )1( 21 1 aa aa M n OLL nnnnnn nn nn bxaxaxa bxaxaxa bxaxaxa =+++ =+++ =+++ L LL L L 2211 22222121 11212111 高斯消元法的第一次消元 )2()2( 2 )2( 2 )2( 2 )2( 22 )2( 22 )1( 1 )1( 12 )1( 121 )1( 11 nnnnn nn nn bxaxa bxaxa bxaxaxa =++ =++ =+++ L LL L L bMAxM 11 = 最终消元形式 = ? AxMMM n 121 L bMMM n 121 L ? 直接法 - 高斯消元法的矩阵表示 第二次消元相当于再左乘单位下三角阵 M 2 bMAxM 11 = bMMAxMM 1212 = 121 , n M MM M M ? =L记 为单位下三角阵 MbUx = UMA=记 0 , )( ≠ k kk a U 对角元素 且为上三角阵 )()( )1( 1 )1( ,11 )1( 1,1 )2( 2 )2( 22 )2( 22 )1( 1 )1( 12 )1( 121 )1( 11 n nn n nn n nn n nnn n nn nn nn bxa bxaxa bxaxa bxaxaxa = =+ =++ =+++ ? ? ? ?? ? ?? LL LL LL MbUx 1? = 直接法-矩阵LU 分解 高斯消元法通过左乘 M,使 MA=U M单位下三角阵， U上三角阵 记 L=M -1 ， L为 单位下三角阵 若 A可逆且顺序主子式不为零，则 A可分解为一个 单位下三角阵 L和一个 上三角阵 U的积 A=LU。 这种分解是唯一的，称矩阵 LU分解。 ),,2,1( 0 )( nk a k kk L= ≠ ),1(,0 1 111 nk aa aa DA kkk k L L MM L =≠=的顺序主子式 若 A可逆，则 存在交换阵 P 使 PA=LU L为单位下三角阵， U为上三角阵。 第 i行与第 k行交换 ),1(00 )()( nkiaa k ik k kk L+=≠= ，但必存在消元中会遇到某个 直接法-矩阵LU 分解 不成立可逆，但顺序主子式若 0≠DA 乘以初等交换阵 UMPAUMA =?= P~交换阵（单位阵经若干次行交换） 直接法 - 对称正定矩阵的分解 正定对称矩阵 A 可分解成对角元素为正的 下三角阵 L 与它的转置矩阵之积，即 T LLA = T LDLA =或 其中 L 是单位下三角阵， D 是元素为正的对角阵。 这种分解称三角分解或 Cholesky 分解。 直接法 - 三对角矩阵的 LU分解 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ????? n nn nnn nnn u cu cu cu l l l ba cba cba cb A 11 22 11 3 2 111 222 11 1 1 1 1 OO O OOOO 11 1 1 2,3, , 2,3, , i i i iiii ub a lin u ublc i n ? ? =? ? ? == ? ? =? =? ? L L 在三次样条插值和其它一些计算中，会遇到求解系数矩阵 A具有 三对角形式 的线性方程组，这时 A的 LU分解 （假定分解存在） 可表为 : L和 U的计 算公式为 3 11 1 1 ,2,, / ()/,1,1 iiii nnn iii i yf yfly i n xyu xycxuin ? + = ? ? =? = ? =? ? =? =? ? L L 线性方程组 Ax=f可通过等价的两个三角形 方程组 Ly=f和 Ux=y求解如下 : 直接法 - 三对角矩阵的 LU分解 求解三对角形方程组的追赶法 直接法-误差分析 （ 1,1） 01.201.1 21 =+ xx x 2 2 x 120 201.1 2 21 21 =+ =+ xx xx bAx = ，如果解 x 对 b 或 A 的扰动敏感，就称方程组是 病态的，也称系数矩阵 A 是病态的。 bxA = ? ? ? ? ? ? = 0 2 x ? ? ? ? ? ? = 1 1 x ? ? ? ? ? ? = 01.2 2 b ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? 2 2 01.1 1 1 1 2 1 x x x对 b的扰动 敏感 (2,0) 向量和矩阵的范数 度量向量、矩阵大小的数量指标 向量范数 最常用的向量范数是 2-范数 2/122 1 2 )( n xxx ++= L xxxx T n ，范数记作设 ),( 1 L= n xxx ++=? L 1 1 1 范数 矩阵范数 AaA nnij ，范数记作设 × = )( 范数 ?= 2)(|||| max2 AAA T λ max λ 表示最大特征根 范数）?= ∑ = 1( max|||| 1 1 n i ij j ||aA 范数）?∞= ∑ = ∞ ( max|||| 1 n j ij i ||aA ),max( 1 n xxx L=?∞ ∞ 范数 向量和矩阵范数的相容性条件 xAAx ?≤ bbxxA δδ +=+ )( bxA δδ = xAb ?≤ b b AA x x δδ ??≤ ?1 1）设 b有扰动 δb，分析 x的误差 δx A的条件数越大,(由 b的扰动引起的) x的误差可能越大 条件数与误差分析 bAx = xAAx ?≤ bAx δδ 1? = bAx δδ ?≤ ?1 Abx /≥ 1 ()A Cond A A A ? =?的 为条定义 件数 b b Acond δ )(= bxxAA =++ ))(( δδ x的 (相对 )误差不超过 b的 (相对 )误差的 Cond(A)倍 , 也大致上是 A的 (相对 )误差的 Cond(A)倍。 条件数与误差分析 2）设 A有扰动 δA，分析 x的误差 δx bAx = A的条件数越大,(由 A的扰动引起的) x的误差越大 A A Acond A A AA x x δδδ )( 1 =??? ? 条件数大的矩阵是 病态矩阵 迭代法---一个例子 123 123 12 3 10 3 14 210 3 5 31014 xxx xxx xx x + += ? +=? ++ = 4.13.01.0 5.03.02.0 4.11.03.0 )( 2 )( 1 )1( 3 )( 3 )( 1 )1( 2 )( 3 )( 2 )1( 1 +??= ++= +??= + + + kkk kkk kkk xxx xxx xxx 0 )0( 3 )0( 2 )0( 1 === xxx L,2,1,0=k 9906.0,9645.0,9906.0 )4( 3 )4( 2 )4( 1 === xxx 4.13.01.0 5.03.02.0 4.11.03.0 213 312 321 +??= ++= +??= xxx xxx xxx 4.1,5.0,4.1 )1( 3 )1( 2 )1( 1 === xxx 1 321 === xxx精确解 4 迭代法 - 雅可比（Jacobi）迭代 将 A 分解为 ULDA ??= ，其中 ),,( 2211 nn aaadiagD L= ， ? ? ? ? ? ? ? ? ? ? ? ? ?= ? ? ? ? ? ? ? ? ? ? ? ? ?= ? ? 0 0 0 , 0 0 0 1,1 112 1,21 21 n n nnnn a aa U aaa a L O MO L L OOM 设对角阵 D非奇异（即 ),1,0 nia ii L=≠ bDf ULDB 1 1 1 1 )( ? ? = +=记 )2,1,0( 1 )( 1 )1( L=+= + kfxBx kk 迭代格式 bAx = bxULDx =+? )( bDxULDx 11 )( ?? ++= 迭代法 - 高斯-塞德尔(Gauss-Sedeil)迭代 在 D非奇异的假设下 )( LD? 可逆，于是得到 bLDfULDB 1 2 1 2 )(,)( ?? ?=?= )2,1,0( 2 )( 2 )1( L=+= + kfxBx kk Jacobi迭代公式 bUxLxDx kkk ++= + )()()1( bUxLxDx kkk ++= ++ )()1()1( Gauss-Seideil迭代公式 4.13.01.0 5.03.02.0 4.11.03.0 )( 2 )( 1 )1( 3 )( 3 )( 1 )1( 2 )( 3 )( 2 )1( 1 +??= ++= +??= + + + kkk kkk kkk xxx xxx xxx 改进 4.13.01.0 5.03.02.0 4.11.03.0 )1( 2 )1( 1 )1( 3 )( 3 )1( 1 )1( 2 )( 3 )( 2 )1( 1 +??= ++= +??= +++ ++ + kkk kkk kkk xxx xxx xxx 迭代法的收敛性 原方程组的解 * x 满足： fBxx += ** 1 )( 1 )1( fxBx kk += + Jacobi迭代 2 )( 2 )1( fxBx kk += + Gauss-Seideil迭代 bDf ULDB 1 1 1 1 )( ? ? = += bLDf ULDB 1 2 1 2 )( )( ? ? ?= ?= fBxx kk += + )()1( 一般迭代形式 )( *)0(*)( xxBxxk kk ?=?次得到迭代 () * () k xxk→→∞序列收敛 的 充要条件 )(0 ∞→→ kB k ?B的所有特征根（取模）小于 1 的特征根是（ 的谱半径 Bni BB i i ni ),1 max)( 1 L= = ≤≤ λ λρ 1)( <Bρ 2）若 A 对称正定，则高斯—塞德尔迭代收敛; 迭代法的收敛性 ||||)( BB ≤ρ谱半径性质： 其中 ||B||是任何一种矩阵范数 () * () k xxk→→∞序列收敛 的 充分条件 赛德尔迭代均收敛；则雅可比和高斯 是严格对角占优的，即若 ? => ∑ ≠ ),,1()1 nibbA ij ijii L 收敛则迭代公式 fBxx kk += + )()1( ，若 1)3 <= qB 越小收敛越快，且 qxx q q xx kkk )()1(*)1( 1 ? ? ≤? ++ 迭代法 -超松弛 SOR 迭代 Gauss-Seideil迭代公式 (1) 1 (1) () 1 (1) () kkk k xDLxUxDbx +?+ ? + =++%null )()1()1( )1( ~ kkk xxx ωω ?+= ++ 改进 (1)k x ω + % k 和x 加权 作平均 1ω > 1ω < 1ω = 超松弛迭代 低松弛迭代 Gauss-Seideil迭代 (1) () 1 1 , ()[(1)], () kk xBxf BDLU D fDLb ωω ω ω ωω ω ωω + ? ? =+ =? +? =? 迭代法 -超松弛 SOR 迭代 收敛充要条件若 A对称正定 02ω< < SOR 迭代 ------解大型稀疏矩阵方程组 5 超定线性代数方程组的最小二乘解 实验 1§ 2.1汽车刹车距离 建立了刹车距离 d与车速 v 的关系： 2 21 vkvkd += 方程个数超过了未知数个数，称为 超定方程组 nivkvkd iii ,,2,1, 2 21 L=+= 数据拟合 已知一组数据，即平面上 n个点( x i ,y i ), i=1,…n, 寻求一个函数 （曲线） y=f(x), 使 f(x)在某种准则下与所有数据点最为接近， 即曲线拟合得最好。 数据拟合问题的提法 (#) + + + + + + + + + x y y=f(x) (x i ,y i ) δ i 使点( x i ,y i ) 与曲线 y=f(x) 的距离 δ i 尽量小 ,i=1,…n 曲线拟合 最小二乘准则 δ L ii i 使f(x )与y(i=1,2, ,n)之差的平方和 (即图中 的平方和)最小 先选定一组函数 r 1 (x), r 2 (x), …r m (x), m<n, 令 f(x)=a 1 r 1 (x)+a 2 r 2 (x)+ …+a m r m (x) （ 1） 其中 a 1 ,a 2 , …a m 为待定系数。 曲线拟合 --线性最小二乘法 的基本思路 记 22 12 11 2 11 (, , ) [() ] [()](2) nn mi ii ii nm kk i i ik Ja a a f x y ar x y δ == == == ? =? ∑∑ ∑∑ L 按照 最小二乘准则 , 问题归结为：求 a 1 ,a 2 , …a m 使 J(a 1 ,a 2 , …a m ) 最小。 线性最小二乘法的求解 ),1( 0 mk a J k L= = ? ? )3( 0])()[( 0])()[( 11 11 1 ? ? ? ? ? ? ? ? ? =? =? ? ∑∑ ∑∑ == == n i m k iikkim n i m k iikki yxraxr yxraxr LL ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? = nmnmn m y y y a a a xrxr xrxr R MM L LL L 11 1 111 ,, )()( )()( 记 )4()()3( yRaRR TT =? )2(])([),,( 2 11 21 iik n i m k km yxraaaaJ ?= ∑∑ == L )4()( yRaRR TT = 当 R T R 可逆时(4)有唯一解 )5()( 1 yRRRa TT ? = 线性最小二乘法的求解 mn nmn m xrxr xrxr R × ? ? ? ? ? ? ? ? ? ? = )()( )()( 1 111 L LL L 问题 1.怎样选择 { r 1 (x), …r m (x)}，以保证系数 {a 1 ,…a m }有唯一解？ {a 1 ,…a m }有唯一解 ← R T R可逆 ←Rank(R T R)=m ← Rank(R)=m ← R列满秩 ← {r 1 (x), …r m (x)}在 x i 点线性无关 (i=1,…n) 分 析 2.为什么要规定 m<n?若 m=n或 m>n,会如何? 若 m>n， a有无穷多解 若 m<n， 超定方程组 , a无解；求 最小二乘解 若 m=n， R可逆时 a有唯一解 强令f(x i )= a 1 r 1 (x i )+ …+a m r m (x i )= y i (i=1, …n) (即曲线 f(x)过全部数据点，此时 J=0) 得 mn nmn m xrxr xrxr R × ? ? ? ? ? ? ? ? ? ? = )()( )()( 1 111 L LL L ,yRa= T n T m yyy aaa ],[ ],[ 1 1 L L = = 问题 分 析 3. 线性最小二乘中的“线性”指的是什么 问题 f(x)对 a线性 ，于是求解 线性 方程组 yRaRR TT =)( 6 线性最小二乘拟合中函数 {r 1 (x), …r m (x)}的选取 1. 通过机理分析建立数学模型来确定 f(x) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + f=a 1 +a 2 x 2. 将数据 (x i ,y i ) i=1, …n 作图，通过直观判断确定 f(x) f=a 1 +a 2 x+a 3 x 2 f=a 1 +a 2 x+a 3 x 2 f=a 1 +a 2 /x f=a 1 exp(a 2 x) f=a 1 exp(a 2 x) f 对 a 非线性，怎么办 1. 求解 bAx= 用左除： bAx \= 。 [x,y]=lu(A) 若 A可逆且顺序主子式不为零， 输出 x 为 单位下三角阵 L，y 为上三角阵U，使 LUA= ； 若A 可逆， x 为一交换阵与单位下三角阵之积 . 2. 矩阵 LU分解 线性方程组数值解法的MATLAB实现 若 A为可逆方阵，输出原方程的解 x 若 A为 n×m矩阵（ n>m），且 A T A可逆，输出原 方程的最小二乘解 x u =chol(A) 对正定对称矩阵A的Cholesky 分解， 输出u 为上三角阵 U，使 A=U T U [x,y,p]=lu(A) 若 A 可逆，输出 x为单位 下三角阵 L，y为上三角阵 U，p 为一交换阵 P， 使 LUPA= . 线性方程组数值解法的MATLAB实现 123 123 12 3 10 3 14 210 3 5 31014 xxx xxx xx x + += ? +=? ++ = 例 . 解 A=[10 3 1;2 -10 3;1 3 10], b=[14 -5 14]', 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分解 MATLAB 5.3.lnk shiyan51 若第 1个方程改为 3x 2 +x 3 =14 结果如何 4． Hilbert 矩阵 ： h=hilb(n) 输出 h为 n 阶 Hilbert 矩阵 3. 范数 条件数 特征值 n=norm(x) 输入 x 为向量或矩阵，输出为 x 的 2-范数 c=cond(x) 输入 x 为矩阵 , 输出为 x 的 2-条件数 r=rcond(x) 输入 x 为方阵 , 输出为 x 条件数倒数 e=eig(x) 输入 x 为矩阵，输出 x 的全部特征值 ? ? ? ? ? ? ? ? ? ? ? ? ?+ + = )12/(1)1/(1/1 )1/(13/12/1 /12/11 nnn n n H L LL LL LL 当 n很大时 Hilbert 矩阵呈病态 线性方程组数值解法的MATLAB实现 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] T MATLAB 5.3.lnk shiyan52 x x1 1.0e+003 * 0.0050 0.0680 -0.1200 -1.3800 0.6300 6.3000 -1.1200 -9.9400 0.6300 5.0400 cond(H)=4.7661e+005 7 y=triu(x) 输入矩阵 x，输出 v 是 x 的上三角阵； v=tril(x) 输入矩阵 x，输出 v 是 x 的下三角阵； 1. 提取（产生）对角阵 v=diag(x) 输入向量x，输出v是以x为对角元素的对 角阵；输入矩阵x，输出v是x的对角元素构成的向量； v=diag(diag(x)) 输入矩阵x，输出v是x的对角元 素构成的对角阵，可用于迭代法中从A中提取D。 2. 提取（产生）上（下）三角阵 线性方程组数值解法的MATLAB实现 v=triu(x,1) 输入矩阵 x，输出 v 是 x 的上三角阵， 但对角元素为 0，可用于迭代法中从 A 中提取 U; v=tril(x,-1) 输入矩阵 x，输出 v 是 x 的下三角阵， 但对角元素为 0，可用于迭代法中从 A 中提取 L 。 x T (雅可比 )x T (高斯—塞德尔 ) 0 (0, 0, 0) (0, 0, 0) 1 (1.4, 0.5, 1.4) (1.4, 0.78, 1.026) 2 (1.11, 1.20, 1.11) (1.0634, 1.0205, 0.9875) 3 (0.929, 1.055, 0.929) (0.9951, 0.9953, 1.0019) 4 (0.9906, 0.9645, 0.9906) (1.0012, 1.0008, 0.9996) 123 123 12 3 10 3 14 2103 5 31014 xxx xxx xx x + += ? +=? ++ = 例 . 用迭代法 解 MATLAB 5.3.lnk shiyan53 稀疏矩阵的处理 ~ MATLAB进行大规模计算的优点 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) 8 aa= 0 0 0 0 (2,3) 8 0 8 8 0 输出 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=toc aa=full(a); tic;xx=aa\b;t2=toc y=sum(x) yy=sum(xx) 例 . 分别用稀疏矩阵和满矩阵求解 Ax=b, 比较计算时间 nn A × ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = 41 1 41 141 14 OO O [ ] T nb L,2,1= 设 0 0 t1, t2相差巨大，说明用稀疏矩阵计算的优点 （ y=yy 用于简单地验证两种方法结果的一致） MATLAB 5.3.lnk shiyan54 用 MATLAB作线性最小二乘拟合 1. 作多项式 f(x)=a 1 x m + …+a m x+a m+1 拟合,可利用 已有程序 : a=polyfit(x,y,m) 输入:数据 x,y (同长度数组); m (拟合多项式次数) 输出:系数 a=[a 1 , …a m ,a m+1 ] (数组)。 2. 对超定方程组 )( 11 nmyaR nmmn <= ××× 仍用 yRa \= 可得最小二乘意义下的解 多项式 在 x点的值 : y=polyval(a,x) 表 1 国民经济各个部门间的关系 分配去向 投入来源 农业 制造业 服务业 外部需求 总产出 农业 15 20 30 35 100 制造业 30 10 45 115 200 服务业 20 60 / 70 150 初始投入 35 110 75 总投入 100 200 150 实例1 投入产出模型 分配去向 投入来源 农业 制造业 服务业 农业 0.15 0.10 0.20 制造业 0.30 0.05 0.30 服务业 0.20 0.30 0 表 2 投入产出表 假定每个部门 的产出与各部 门对它的投入 成正比，得到 投入系数 。 8 4）如果对于任意给定的、非负的外部需求，都能 得到非负的总产出，模型就称为可行的。问为使模 型可行，投入系数应满足什么条件？ 1）设有n个部门，已知投入系数，给定外部需求， 建立求解各部门总产出的模型。 2）设投入系数如表2所给，如果今年对农业、制造 业和服务业的外部需求分别为50，150，100亿元， 问这三个部门的总产出分别应为多少。 3）如果三个部门的外部需求分别增加1个 单位，它们的总产出应分别增加多少。 实例1 投入产出模型 1）基本模型 投入系数矩阵 nnij aA × = )( 产出向量 T n xxx ),( 1 L= 需求向量 T n ddd ),( 1 L= ),,2,1( 1 nixdx i n j iij L==+ ∑ = 平衡关系 x i : 第 i个部门的产出， x ij : 第 i个部门对第 j个部 门的投入， d i : 第 i个部门的外部需求 jijij xxa /= 投入系数 产出 投入 部门 1 部门 i 部门 n 外部需求 总产 出 部门 1 x11 x1i x1n d1 x1 部门 i xi1 xii xin di xi 部门 n xn1 xni xxn dn xn 初始投入 x01 x0i x0n 总投入 x1 xi xn ),,2,1( 1 nixdxa i n j ijij L==+ ∑ = dAxx += dxAI =? )( dAIx 1 )( ? ?= dAIx 1 )( ? ?=基本模型 2）设农业、制造业和服务业的外部需求分别为 50，150，100亿元，求三个部门的总产出。 ? ? ? ? ? ? ? ? ? ? = 03.02.0 3.005.03.0 2.01.015.0 A T d ]10015050[= x=(139.2801,267.6056,208.1377) T MATLAB 5.3.lnk shiyan55 若 ?d=(1,0,0) T , 即农业外部需求增加1单位时，三部门 总产出应分别增加1.3459，0.5634，0.4382单位。 即 C的第1列 。 C的第2，3列给出了什么？ C=1.3459 0.2504 0.3443 0.5634 1.2676 0.4930 0.4382 0.4304 1.2167 3）若三部门的外部需求分别增加1个单位， 求它们的总产出的增量。 dAIx 1 )( ? ?=基本模型 1 )( ? ?= AIC记 当需求增加 ?d时，总产出增量 dCx ?=? 4）如果对于任意给定的、非负的外部需求，都能 得到非负的总产出，模型就称为可行的。问为使模 型可行，投入系数应满足什么条件？ 模型可行 1 ))(( + ?=+++? kk AIAAIAI L 0,)( 1 ≥?= ? AdAIx 其中基本模型 00 ≥?≥? xd 0)( 1 ≥? ? AI )(0 ∞→→ kA k ∑ ∞ = ? =? 0 1 )( k k AAI 是否成立？问 )(0 ∞→→ kA k 0≥A ),1(11 11 1 njxxaA j n i ij n i ij L=<?<?< ∑∑ == 001 1 →?→?< k k AAA 模型可行 00 ≥?≥? xd BAAB ?≤矩阵范数定义 kk AA ≤ )(0 ∞→→ kA k 成立),1( 1 njxx j n i ij L=< ∑ = 若A是由实际数据得 到，只要初始投入非负 模型可行 jijij xxa /= max|||| 1 1 ∑ = = n i ij j ||aA 9 实验 2给出的模型为 实例2 一年生植物的繁殖 (实验 2) nkqxpxx kkk ,,3,2,0 21 L==++ ?? 121 ,(1)p abc q ab a bc=? =? ? Ax=B ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = pq pq pq pq p A 1 1 1 1 OOO ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? 1 2 3 2 1 n n x x x x x x M 0 0 0 0 n qx B x ??? ?? ?? ?? = ?? ?? ?? ?? ? ?? M 求解 Ax=B可得第二年（及以后诸年）植物的数量 开始有 100棵植物，要求 50年后有 1000棵植物 实例2 一年生植物的繁殖 12 0.5, 0.25, 0.2, 10aa bc= ===设 p=-1;q=-0.05;x0=100;xn=1000; n=49; A1=sparse(1:n,1:n,p,n,n); A2=sparse(1:n-1,2:n,1,n,n); A3=sparse(2:n,1:n-1,q,n,n); A=A1+A2+A3; i=[1,n];j=[1,1];s=[-q*x0,-xn]; B=sparse(i,j,s,n,1); x=A\B; x1=x(1), % 输出第 2年植物数量 101.7097 k=0:n+1;xx=[x0,x',xn]; plot(k,xx),grid, 0 10 20 30 40 50 0 200 400 600 800 1000 k x k 图 4 一年生植物繁殖 x k MATLAB 5.3.lnk shiyan56 实例3 汽车刹车距离 (续 ) 实验 1§ 2.1建立了刹车距离 d与车速 v的关系： 2 21 vkvkd += nivkvkd iii ,,2,1, 2 21 L=+= 数据拟合 (#) 车速与刹车距离的实际数据记作 (vi, di), i=1,2,… ,7（见实验 1表 1) (#)是超定方程组 T n T n nn ddykk vv vv ),(,),(, 121 2 2 2 11 LL == ? ? ? ? ? ? ? ? ? ? =Φ × β yβΦ= 12 0.6522, 0.0853kk= = 结果 MATLAB 5.3.lnk shiyan57 布置实验 目的 1)学会用 MATLAB软件数值求解线性代数方程组，对迭 代法的收敛性和解的稳定性作初步分析； 2)通过实例学习用线性代数方程组解决简化的实际问题。 内容 课上布置，或参见网络学堂