为什么县学习能性才覆 大学数学实验 的教值解法 Experiments in Mathematics 许多实际问题归结为线性(代数)方程组 机械设备、土建结构的受力分析 经济计划 实验5线性代数方程组的数值解法 输电网络、管道系统的参数计算企业管理 大型的方程组需要有效的数值解法 半大歇教科系 教值解法的稳定性和收敛性问题需要注意 实验5的走內内率 线性方程组的一般形式、两类解法 anx1+a2x2+…+ax=b 1.两类数值解法: a2x+a2x2+…+a2x,=b2 或AX=b 直接方法;选代方法 ax.+ax.+…+ax=b 2.超定线性方程组的最小二乘解 经过有限次算术迳算求出精确解(实际上 3.线性方程组数值解法的 MATLAB实现 由于有舍入误差只能得到近似解) 高 (GBu8)消元法及与它密切相关的矩阵UU分解 4.实际问题中方程组的数值解 选代法从初始解出发,根据设计好的步骤用次 求出的近似解逼近精确解 雅可比(J8cobi 选代法和高斯一塞德尔( Gauss-Seide1)选代法 (学静学实鉴 学学实纷 直接法一高斯消元法 直接法-列主元素消元法 a1x1+a1x,+…+anx 高斯消元法条件 a2x2+…+a2xn=2 ak(绝对值)很小时, +ax,= 过 用它作除数会导致舍入误 ami. 差的很大增加 qx1+ax2+……+amxn=b +……+ak)x=kk) 解￥ amx,t++. 办法 条件 过 最大的一个(列主元) 程 将列主元所在行与第k行交换后,再按上面的高斯消元 (k=1,2,…,n) 进行下去,称为列主元素消元法
1 大学数学实验 Experiments in Mathematics 实验5 线性代数方程组的数值解法 清华大学数学科学系 • 许多实际问题归结为线性（代数）方程组 • 大型的方程组需要有效的数值解法 • 数值解法的稳定性和收敛性问题需要注意 为什么要学习线性方程组 的数值解法 机械设备、土建结构的受力分析 输电网络、管道系统的参数计算 经济计划 企业管理 3. 线性方程组数值解法的MATLAB实现 4. 实际问题中方程组的数值解。 1. 两类数值解法: 直接方法；迭代方法 实验5的主要内容 2. 超定线性方程组的最小二乘解 线性方程组的一般形式、两类解法 直接法 经过有限次算术运算求出精确解（实际上 由 于 有 舍 入 误 差 只 能 得 到 近 似 解 ） ---- 高 斯 （Gauss）消元法及与它密切相关的矩阵LU分解 迭代法 从初始解出发，根据设计好的步骤用逐次 求出的近似解逼近精确解 ---- 雅可比（Jacobi） 迭代法和高斯—塞德尔（Gauss—Seidel）迭代法 n n nn n n n n n n a x a x a x b a x a x a x b a x a x a x b + + + = + + + = + + + = L LL L L 1 1 2 2 21 1 22 2 2 2 11 1 12 2 1 1 或 AX=b 直接法---高斯消元法 (2) (2) 2 (2) 2 (2) 2 (2) 2 2 (2) 22 (1) 1 (1) 2 1 (1) 1 12 (1) 11 n nn n n n n n n a x a x b a x a x b a x a x a x b + + = + + = + + + = L LL L L ( ) ( ) ( 1) 1 ( 1) 1 1, ( 1) 1, 1 (2) 2 (2) 2 2 (2) 22 (1) 1 (1) 2 1 (1) 1 12 (1) 11 n n n n nn n n n n n n n n n n n n n n a x b a x a x b a x a x b a x a x a x b = + = + + = + + + = − − − − − − − − LL LL LL 消 元 过 程 回 代 过 程 ( 1,2, , ) 0 ( ) k n a k kk = L ≠ 条件 直接法-列主元素消元法 0( 1,2, , ) ( ) a k n k 高斯消元法条件 kk ≠ = L 解决 办法 ( , ) ( ) a i k n k 选 ik = L 最大的一个（列主元） 将列主元所在行与第k行交换后, 再按上面的高斯消元 法进行下去，称为列主元素消元法。 (绝对值 )很小时， 用它作除数会导致舍入误 差的很大增加 (k) kk a ( ) ( ) ( ) ( ) ( ) ( ) k n n k k nn k nk k n k k k kn k kk a x a x b a x a x b + + = + + = LL LL LL
直接法一高斯消元法的矩阵表示 直接法一高斯消元法的矩阵表示 高斯消元法的第一次消元 第二次消元相当于再左乘单位下三角阵M a11+a12x2+……+a1nx ax+ax2+…+amxn= M4x=Mb口M2MAx=M2Mb a,x,+anx+-+a, x=b a2x2+…+a23xn=b2 最终消元形式M。…M2M1Ax=Mn…MMb 记Mn1…M2M1 an1x1+a2x2+…+amxn=b anx2+…+axn=b M为单位下三角阵 相当于方程 -a/a1 Ux= Mb AX=b两边 左乘单位下 M, Ax= M, b U为上三角阵,且 ap-L- r+ap=, =b el 三角阵M1 对角元素a()≠0 x=U-Mb 直油-矩I分解 直油-矩陈分解 若A可逆,但顺序主子式D≠0不成立 0已A的顺序主子式D (k=1,2,…,n) 消元中会遇到某个an=0,但必存在aA)≠0=k+1…n) 高斯消元法通过左乘M使MA=U记L=M,L为 第行与第A行交换日乘以初等交换阵 M单位下三角阵,U上三角阵 单位下三角阵 若A可逆且顺序主子式不为零,则A可分解为一个 P~交换阵(单位阵经若干次行交换) 单位下三角阵和一个上三角阵U的积A=LU 若A可逆,则存在交换阵P使PA=LU 这种分解是唯一的,称矩阵LU分解。 L为单位下三角阵,U为上三角阵 学学实 (大学数学实验) 直接法一对称正定矩阵的分解 直接法一三对角矩阵的L解 正定对称矩阵A可分解成对角元素为正的 在三次样条插值和其它一些计算中,会遇到求解系数矩阵A具有 下三角阵L与它的转置矩阵之积,即 对角形式的线性方程组,这时A的LU分解(假定分解存在) 可表为: 或A=LDL 4,h c 其中L是单位下三角阵,D是元素为正的对角阵 这种分解称三角分解或 Cholesky分解 L和U的计 算公式为 =b-l-1t=2,3…n
2 直接法 - 高斯消元法的矩阵表示 相当于方程 AX=b 两 边 左乘单位下 三角阵M1 ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = / 1 / 1 1 (1) 11 (1) 1 (1) 11 (1) 21 1 a a a a M n LL O n n nn n n n n n n a x a x a x b a x a x a x b a x a x a x b + + + = + + + = + + + = L LL L L 1 1 2 2 21 1 22 2 2 2 11 1 12 2 1 1 高斯消元法的第一次消元 (2) (2) 2 (2) 2 (2) 2 (2) 2 2 (2) 22 (1) 1 (1) 2 1 (1) 1 12 (1) 11 n nn n n n n n n a x a x b a x a x b a x a x a x b + + = + + = + + + = L LL L L M1Ax = M1b 最终消元形式 M n−1LM 2M1Ax = M M M b n−1L 2 1 直接法 - 高斯消元法的矩阵表示 第二次消元相当于再左乘单位下三角阵M2 M1Ax = M1b M2M1Ax = M2M1b 1 21 , Mn MM M M 记 − L = 为单位下三角阵 记 MA =U Ux = Mb 0 , ( ) ≠ k kk a U 对角元素 为上三角阵 且 ( ) ( ) ( 1) 1 ( 1) 1 1, ( 1) 1, 1 (2) 2 (2) 2 2 (2) 22 (1) 1 (1) 2 1 (1) 1 12 (1) 11 n n n n nn n n n n n n n n n n n n n n a x b a x a x b a x a x b a x a x a x b = + = + + = + + + = − − − − − − − − LL LL LL x U Mb −1 = 直 接 法 - 矩 阵 LU 分 解 高斯消元法通过左乘M,使MA=U M单位下三角阵，U上三角阵 记 L=M-1，L为 单位下三角阵 若A可逆且顺序主子式不为零，则A可分解为一个 单位下三角阵L和一个上三角阵U的积 A=LU。 这种分解是唯一的，称 矩阵LU分解。 ( 1,2, , ) 0 ( ) k n a k kk = L ≠ 0, ( 1, ) 1 11 1 k n a a a a A D k kk k L L M M L 的顺序主子式 = ≠ = 若A可逆，则存在交换阵 P 使 PA=LU L为单位下三角阵，U为上三角阵。 第i行与第k行交换 0 0( 1, ) ( ) ( ) a a i k n k ik k 消元中会遇到某个 kk = ，但必存在 ≠ = + L 直 接 法 - 矩 阵 LU 分 解 若A可逆，但顺序主子式 D ≠ 0不成立 乘以初等交换阵 MA =U ⇒ MPA =U P~交换阵（单位阵经若干次行交换） 直接法 - 对称正定矩阵的分解 正定对称矩阵 A 可分解成对角元素为正的 下三角阵 L 与它的转置矩阵之积，即 T A = LL T 或 A = LDL 其中 L 是单位下三角阵，D 是元素为正的对角阵。 这种分解称三角分解或 Cholesky 分解。 直接法 - 三对角矩阵的LU分解 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = − − − − − n n n n n n n n n u u c u c u c l l l a b a b c a b c b c A 1 1 2 2 1 1 3 2 1 1 1 2 2 2 1 1 1 1 1 1 O O O O O O O 1 1 1 1 2,3, , 2,3, , i i i i i ii u b a l in u u b lc i n − − ⎧ = ⎪ ⎪ ⎨ = = ⎪ ⎪⎩ =− = L L 在三次样条插值和其它一些计算中，会遇到求解系数矩阵A具有 三对角形式的线性方程组，这时A的LU分解（假定分解存在） 可表为: L和U的计 算公式为
A x=b 大学静学实扮 直接法一三对角矩阵的L解 真染-误基分「1「x 线性方程组Ax=通过等价的两个三角形 X+x2=2 11.01x2 方程组Ly=/和Ux=y求解如下 x1+10lx=2 y=f x+1.01x2=2.01 y=f-y,i=2…n x对b的扰动 jx x=(y-cx)l4,i=n-1-…1 Ax=b,如果解x对b或A的扰动敏感,就称方程组是 求解三对角形方程組的追赶法 病态的,也称系数矩阵A是病态的 向量和矩阵的救底量向量、绳阵大小的量指 条件敷与误差分析Ax=bA|s4| 向量范数设x=(x,…xn),范数记作| A(x+x)=b+荡口Aox=b 最常用的向量范数是2-范数2=(x+ 日科H网 1-范数风=x+…+∞-范数|L=max(…|x |4叫|≥|4 矩阵范款设A=(an)nn,范数记作|4 ‖北l= YmaNA)2-范数max表示最大特征根 1Ah=maxa(-范数)AL2=maxl(∞-范数 定义的条件数为Cond(4)=|4-1|4 A的条件数越大,(由的扰动引起的)x的误差可能越大 量和细降花的幅亭|4≤|4 (学静学实鉴 (大学数学实验) 条件敷与误差分析 选代法一一个例子 2)设A有扰动6A,分析x的误差8x x1=-0.3x2-0.x3+1.4 (A+SA(x+ Sx)=b 0.lx1-0.3x+14 cond(a) k=0,1,2 A的条件数越大,(由A的扰动引起的)x的误差越大 x3(k+)=-01x)-03x2()+14 x的(相对)误差不超过b的(相对)误差的ConA)倍, x0)=x0)=0x=1.4,x2=0.5,x 也大致上是A的(相对误差的Con4)倍 x1)=0.9906x=09645x4=0 条件教大的矩阵是病态矩阵 精确解x1=x2=x3=1
3 1 1 1 1 , 2, , / ( )/ , 1, ,1 i i ii n nn i i ii i y f y f ly i n x yu x y cx u i n − + ⎧ = ⎨ ⎩ =− = ⎧ = ⎨ ⎩ = − =− L L 线性方程组Ax=f可通过等价的两个三角形 方程组Ly=f和Ux=y求解如下 : 直接法 - 三对角矩阵的LU分解 求解三对角形方程组的追赶法 直接法-误差分析 （1,1） 1.01 2.01 1 2 x + x = x2 2 x 0 2 1 1.01 2 2 1 2 1 2 + = + = x x x x Ax = b，如果解 x 对 b 或 A 的扰动敏感，就称方程组是 病态的，也称系数矩阵 A 是病态的。 A x = b ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 0 2 x ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 1 1 x ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 2.01 2 b ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 2 2 1 1.01 1 1 2 1 x x x对b的扰动 敏感 (2,0) 向量和矩阵的范数 度量向量、矩阵大小的数量指标 向量范数 最常用的向量范数是 2-范数 2 2 1/ 2 2 1 ( ) n x = x +L+ x x x x x T 设 = ( 1 ,L n ) ，范数记作 n − x = x +L+ x 1 1 1 范数 矩阵范数 设 A = (aij ) n×n，范数记作 A ||A||2= max(A A) 2 −范数 T λ λmax 表示最大特征根 = ∑ −范数） = || || max (1 1 1 n i ij j A |a | = ∑ ∞ −范数） = || ||∞ max ( 1 n j ij i A |a | max( , ) 1 n ∞ − x = x L x 范数 ∞ 向量和矩阵范数的相容性条件 Ax ≤ A ⋅ x A(x + δx) = b + δb Aδx = δb b ≤ A ⋅ x b b A A x δx δ ≤ ⋅ ⋅ −1 1）设b有扰动 δb，分析x的误差 δx A的条件数越大,(由b的扰动引起的)x的误差可能越大 条件数与误差分析 Ax = b Ax ≤ A ⋅ x δx A δb −1 = δx ≤ A ⋅ δb −1 x ≥ b / A 1 A Cond A A A ( ) − 定义 的条件数为 = ⋅ b b cond A δ = ( ) ( A + δA)( x + δx) = b x的(相对)误差不超过b的(相对)误差的Cond(A)倍, 也大致上是A的(相对)误差的Cond(A)倍。 条件数与误差分析 2）设A有扰动 δA，分析x的误差 δx Ax = b A的条件数越大,(由A的扰动引起的)x的误差越大 A A cond A A A A A x δx δ δ ( ) 1 ≅ ⋅ ⋅ = − 条件数大的矩阵是病态矩阵 迭代法 --- 一 个 例 子 1 23 1 23 12 3 10 3 14 2 10 3 5 3 10 14 x xx x xx xx x + + = − + =− ++ = 0.1 0.3 1.4 0.2 0.3 0.5 0.3 0.1 1.4 ( ) 2 ( ) 1 ( 1) 3 ( ) 3 ( ) 1 ( 1) 2 ( ) 3 ( ) 2 ( 1) 1 = − − + = + + = − − + + + + k k k k k k k k k x x x x x x x x x 0 (0) 3 (0) 2 (0) 1 x = x = x = k = 0,1,2,L 0.9906, 0.9645, 0.9906 (4) 3 (4) 2 (4) x1 = x = x = 0.1 0.3 1.4 0.2 0.3 0.5 0.3 0.1 1.4 3 1 2 2 1 3 1 2 3 = − − + = + + = − − + x x x x x x x x x 1.4, 0.5, 1.4 (1) 3 (1) 2 (1) x1 = x = x = 精确解 x1 = x2 = x3 = 1
逸代浩-雅可比( Jacobi)迭代 地代-高斯-塞德尔( Gauss- Seder1)迭代 将A分解为A=D-L-U,其中D=diag(a1 cobi选代公式D(+=Lx(k)+Ux()+b 0 0 )+03x4+05 uA1203“+1进x=02x+40x”+05 设对角阵D非奇异(即a≠0,i=1…n)Ax=b Gauss-Seideily选代公式Dx4+)=Lxk+)+Ux)+b 4 Dx-(L+U)x=b x=D(L+U)x+D-b 在D非奇异的假设下(D-L)可逆,于是得到 记B1=D(L+U B2=(D-L)U, f2=(D-L)b 迭代格式 x(+=B2x()+/2(k=0,2…) f=Dbx+=Bx+f(k=012…) 代沽的收笕 B=D-(L+L 选代法的收做性 Jacobi选代x=Bx)+f1f=Db 序列收敛x*)→x'(k→∞)的充分条件 Gauss-Seideili选代x=B2x1)+f2B=(D-0 1)若A是严格对角占优的,即>∑内(=1…m) 一般选代形式x=Bx4+f f2=(D-L)-b 则雅可比和高斯一赛德尔迭代均收敛 星组的解x满足:x=Bx+f 代k次得到x“-x=B(x-x) 2)若A对称正定,则高斯一塞德尔迭代收敛 序列收敛x→x(k→∞)的充要条件 3)若B=q1 <1 超松弛迭代 低松弛迭代 Gaus- Seidel1遗代 SOR选代—解大型稀疏矩阵方程組
4 迭代法 - 雅可比（Jacobi）迭代 将 A 分解为 A = D − L − U ，其中 ( , , ) D = diag a11 a22 Lann ， ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = − ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = − − − 0 0 0 , 0 0 0 1,1 12 1 1 2 , 1 21 n n n n n n a a a U a a a a L O O M L L M O O 设对角阵D非奇异（即a 0, i 1, n) ii ≠ = L f D b B D L U 1 1 1 1 ( ) − − = 记 = + ( 0,1,2 ) 1 ( ) 1 x (k+1) =Bx k + f k = L 迭代格式 Ax = b Dx − (L +U )x = b x D L U x D b 1 1 ( ) − − = + + 迭代法 - 高斯-塞德尔(Gauss-Sedeil)迭代 在D非奇异的假设下(D − L)可逆，于是得到 B D L U f D L b 1 2 1 2 ( ) , ( ) − − = − = − ( 0,1,2 ) 2 ( ) 2 x(k+1) = B x k + f k = L Jacobi迭代公式 Dx Lx Ux b k k k = + + ( +1) ( ) ( ) Dx Lx Ux b k k k = + + Gauss-Seideil迭代公式 ( +1) ( +1) ( ) 0.1 0.3 1.4 0.2 0.3 0.5 0.3 0.1 1.4 ( ) 2 ( ) 1 ( 1) 3 ( ) 3 ( ) 1 ( 1) 2 ( ) 3 ( ) 2 ( 1) 1 =− − + = + + =− − + + + + k k k k k k k k k x x x x x x x x x 改进 0.1 0.3 1.4 0.2 0.3 0.5 0.3 0.1 1.4 ( 1) 2 ( 1) 1 ( 1) 3 ( ) 3 ( 1) 1 ( 1) 2 ( ) 3 ( ) 2 ( 1) 1 = − − + = + + = − − + + + + + + + k k k k k k k k k x x x x x x x x x 迭代法的收敛性 原方程组的解 *x 满足： x = Bx + f * * 1 ( ) 1 ( 1) x B x f k k = + + Jacobi迭代 2 ( ) 2 ( 1) x B x f k k = + + Gauss-Seideil迭代 f D b B D L U 1 1 1 1 ( ) − − = = + f D L b B D L U 1 2 1 2 ( ) ( ) − − = − = − x Bx f k k = + ( +1) ( ) 一般迭代形式 ( ) ( ) * (0) * k x x B x x k k 迭代 次得到 − = − () * ( ) k 序列收敛 的 x xk → →∞ 充要条件 B → 0(k → ∞) k ⇔B的所有特征根（取模）小于1 （ 是 的特征根 的谱半径 i n B B B i i i n 1, ) ( ) max 1 = L = ≤ ≤ λ ρ λ ρ(B) ∑ = ≠ 1) A b b (i 1, n), j i ii ij L 3) 若 B = q 1 ω < 1 ω = 1 超松弛迭代 低松弛迭代 Gauss-Seideil迭代 ( 1) ( ) 1 1 , ( ) [ (1 ) ], ( ) k k x Bx f B DL U D f D Lb ω ω ω ω ωω ω ω ω + − − = + = − +− = − 迭代法 -超松弛 SOR 迭代 若A对称正定 收敛充要条件 0 2 < ω < SOR 迭代------解大型稀疏矩阵方程组
大学静学实扮 鬼定觉代歙方翟组的录小二最 曲线拟合 使点(xy)与曲线y=/x) 的距高8尽量小户1,n 实验1§2.1汽车刹车距高建立了刹车距离d与车速v f(r) 的关系 d,=k1v,+k2y2,i=1,2,…,n 方程个数超过了未知数个数,称为超定方程组 最小二乘准则 已知一组数据,即平面上n个点x),p=1n寻求一个函数 (曲线)y=fx).使(x)在某种准则下与所有数据点最为接近 使f(x1)与y1(i=1,2…,n)之差的平方和 即曲线拟合得最好 (即图中δ的平方和)最小 曲线拟合-黛恤量小二晁油的基本思路 线性最小二乘法的求解 先选定一组函数rA(x),r2(x),,rm(x),mn,会如何? r1(x1)…m(x1) 当RR可逆时(4)有唯一解R 分强令f(x)=a(x)+…+ am(F(G=1,n a=(RR r1(xn)…r(xn) r1(x1)…rn( a={a1,…an Ra=y y1,…ynj 同题 1.怎样选择Gr(x,…rn(o),以保证系数 {apa)有唯一解? 若m>m,a有无穷多解若m=,R可逆时a有唯一解 {a1,a有唯一解←RR可逆←Rank(RR)=m 若m<n,超定方程组,a无解;求最小二乘解 ←Rank(R)=m←R列满秩 问题3.线性最小二乘中的“线性”指的是什么 ←{r(x),…,rm(x)}在x点线性无关(产=1,…n) f(x)对a线性,于是求解线性方程组(RR=Ry
5 超定线性代数方程组的最小二乘解 实验1§2.1汽车刹车距离建立了刹车距离d与车速v 的关系： 2 1 2 d = k v + k v 方程个数超过了未知数个数，称为超定方程组 d k v k v i n i i i , 1,2, , 2 = 1 + 2 = L 数据拟合 已知一组数据，即平面上 n个点(xi ,yi ), i=1,…n, 寻求一个函数 （曲线）y=f(x), 使f(x)在某种准则下与所有数据点最为接近， 即曲线拟合得最好。 数据拟合问题的提法 (#) + + + + + + + + + x y y=f(x) (xi ,yi ) δi 使点(xi ,yi ) 与曲线 y=f(x) 的距离δi尽量小,i=1,…n 曲线拟合 最小二乘准则 δ i i L i 使 f(x )与 y (i=1,2, ,n)之差的平方和 (即图中 的平方和)最小 先选定一组函数 r1(x), r2(x), …rm(x), mn,会如何? 若m>n，a有无穷多解 若m<n，超定方程组, a无解；求最小二乘解 若m=n，R可逆时a有唯一解 强令f(xi)= a1r1(xi )+ …+amrm(xi )= yi (i=1, …n) (即曲线 f(x)过全部数据点，此时J=0) 得 n m n n m m r x r x r x r x R × ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ( ) ( ) ( ) ( ) 1 1 1 1 L L L L Ra= y, T n T m y y y a a a [ , ] [ , ] 1 1 L L = = 问题 分 析 问题 3. 线性最小二乘中的“线性”指的是什么 f(x)对a线性，于是求解线性方程组 R R a R y T T ( ) =
学静学笑) 大学酸学笑验 线性最小二乘拟合中函数{r/(x),…rm(x)}的选取 线性方程组数值解法的 MATLAB实现 1.通过机理分析建立教学模型来确定几x) 1.求解Ax=b用左除:x=A\b。 将效据(x)P1…n作图,通过直观判断确定jx faar+axr 若A为可逆方阵,输出原方程的解x -ar+ax+ar 若A为mXm矩阵(m>m),且∮A可逆,输出原 方程的最小二乘解 2矩阵LU分解 exp(axt) x,y]=1u(A)若A可逆且顺序主子式不为零,输出x为 单位下三角阵L,y为上三角阵U,使A=LU;若A可逆, x为一交换阵与单位下三角阵之积 ∫对a非线性,怎么办 线性方程组数值解法的 MATLAB实现 10x1+3x2+x3=14 x1+3x2+10x2=14 x,y,p]=1u(A)若A可逆,输出x为单位 下三角阵L,y为上三角阵U,p为一交换阵P 并对系数矩阵=[1031;2-103:1310 使PA=LL 作LU分解 u=ch01)对正定对称矩阵A的les分解, [L1,U1]=1u(A) 着第1个方程改为 A1=L1+U1 输出u为上三角阵U,使A [L2,U2,P]=1u(A) 结果如何 A2=L2+U2 (学静学实鉴 大学酸学实验) 线性方程组数值解法的 MATLAB实现 少观察 Hilbert矩阵的病态性 3.慕数条件微特征值 =norm(x)输入x为向量或矩阵,输出为x的2-范数 例.Hx=b,其中 c=cond(x)输入x为矩阵,输出为x的2-条件数 r= cond(x)输入x为方阵,输出为x条件数倒数 输入x为矩阵,输出x的全部特征值 H=hilb(5) h=rats(H), 4. Hilbert矩阵:h=hilb(n) b=ones(5, 1) 0.00500.0680 输出h为n阶 Hilbert矩阵 -0.1200 b(5)=1.1 1/n x1=H\ 1/21/3 当n很大时 Hilbert x,x1], 0.6300 矩阵呈病态 nl=cond (H) n2=cond (H), cnd(H)=4.7661e+005 l/nl/n+1)…1/O2n-1)
6 线性最小二乘拟合中函数 {r1(x), …rm(x)}的选取 1. 通过机理分析建立数学模型来确定 f(x) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + f=a1+a2x 2. 将数据 (xi ,yi ) i=1, …n 作图，通过直观判断确定 f(x) f=a1+a2x+a3x2 f=a1+a2x+a3x2 f=a1+a2/x f=a1exp(a2x) f=a1exp(a2x) f 对a 非线性，怎么办 1. 求解Ax=b 用左除： x = A\b 。 [x,y]=lu(A) 若 A 可逆且顺序主子式不为零， 输出 x 为 单位下三角阵 L，y 为上三角阵 U，使 A=LU； 若 A 可逆， x 为一交换阵与单位下三角阵之积. 2. 矩阵LU分解 线性方程组数值解法的MATLAB实现 若A为可逆方阵，输出原方程的解x 若A为n×m矩阵（n>m），且ATA可逆，输出原 方程的最小二乘解x u =chol(A) 对正定对称矩阵 A 的 Cholesky 分解， 输出 u 为上三角阵 U，使 A=UT U [x,y,p]=lu(A) 若 A 可逆，输出 x 为单位 下三角阵 L，y 为上三角阵 U，p 为一交换阵 P， 使PA=LU. 线性方程组数值解法的MATLAB实现 1 23 1 23 12 3 10 3 14 2 10 3 5 3 10 14 x xx x xx 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个方程改为 3x2+x3=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 的全部特征值 ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + − + = 1/ 1/( 1) 1/(2 1) 1/ 2 1/3 1/( 1) 1 1/ 2 1/ n n n 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
少线性方程组数值解法的 MATLAB实现 v=triu(x,1)輸入矩阵x,输出v是x的上三角阵 1.提取(产生)对角阵 对角元永为0,可用于选代法中从A中提取歌 v=diag(x)入向量x,输出v是以x为对角元素的对 v=tri1(x,-1)輸入矩阵x,输出v是x的下三角阵 角阵;输入矩阵x,輸出v是x的对角元素构成的向量; 对角元素为0,可用于选代法中从A中提取L 例.用迭代法 +3+= v=diag(diag(x))輸入矩阵x,输出v是x的对角元 2-+x=5 shiyan 几起53】nh 构成的对角阵,可用于选代法中从A中提取D x+3+=14 2.提取(产生)上(下)三角阵 雅可比 高斯一嘉德尔 y=triu(x)輸入矩阵x,輸出V是x的上三角席; 0,1,11) v=tril(x) 输入矩阵x,输出V是x的下三角阵 (09906,0.9645,0.9906) 种疏矩阵的处理~ MATLAB进行大规模计算的优点 例.分别用稀疏矩阵和满矩阵求解Ac=b,比较计算时间 a= sparse(x,c,v,m,n)在第x行、第c列输 =500;b=[1:n] 入数值v,矩阵共m行n列,输出a为稀疏矩阵,只 141 0 1l=sparse(l: n, 1 4,n,n) 给出(x,∞)及v a2=sparse(2: n, 1: n-1, 1, n, n) tic: x=a\b: tl=toc aa=fu11(a)输入种疏矩阵a,输出aa为满矩阵 a=full(a) (包含零元) tic: xx=aa\b: t2=toc sum(x) a= sparse(2,2:3,8,2,4),aa=fu11(a), 2 hivan5s-4 2,3)8 t,t相差巨大,说明用稀疏矩阵计算的优点 (y=yy用于简单地验证两种方法结果的一致) (学静学实鉴 (大学数学实验) F用 MATLAB1作线性最小二乘拟合 实例1投入产出模型 国民经济名个邮门间的关系 1.作多项式∫x)=a1m+…+anx+anm拟合可利用 造业最务业外部霜求总产 已有程序:a= polyfit(x,y,m 输入:数据x,y(同长度数组);m(拟合多项式次数) 输出:系数a=[a1,-an,a+1(数组)。 多项式在x点的值:y= polya1(a,x) 假定每个部门 2.对超定方程组 R a=yna1(m<n 的产出与各部 农业新造业最务业 门对它的投入 仍用a=R\y可得最小二乘意义下的解 成正比,得到 投入系数
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。 xT(雅可比) xT(高斯—塞德尔) 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) 1 23 1 23 12 3 10 3 14 2 10 3 5 3 10 14 x xx x xx 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, 比较计算时间 n n A × ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 1 4 1 1 4 1 4 1 4 1 O O O [ ] T b = 1,2 ,L n 设 0 0 t1, t2相差巨大，说明用稀疏矩阵计算的优点 （y=yy 用于简单地验证两种方法结果的一致） MATLAB 5.3.lnk shiyan54 用MATLAB作线性最小二乘拟合 1. 作多项式 f(x)=a1xm+ …+amx+am+1拟合,可利用 已有程序: a=polyfit(x,y,m) 输入:数据x,y (同长度数组);m (拟合多项式次数) 输出:系数a=[a1, …am , am+1] (数组)。 2. 对超定方程组 ( ) Rn×mam×1 = yn×1 m < n 仍用 a = R \ y 可得最小二乘意义下的解 多项式在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 投入产出表 假定每个部门 的产出与各部 门对它的投入 成正比，得到 投入系数
物1入产出桃型 X学酸学实 1)基本模型 )设有n个鄣门,已知投入系敷,给定外部需求 x第个部门的产出 建立求解各部门总产出的模型。 x第个部门对第f个部 2)设投入系数如表2所给,如果今年对农业、制造 d;第个邮门的外部需求 业和服务业的外部需求分别为50,150,100亿元 问这三个部门的总产出分别应为多少。 平衡关系 3)如果三个部门的外鄣需求分别增加1个 xy +a== 口2 单位,它们的总产出应分别增加多少 4)如果对于任意给定的、非负的外部需求,都能 投入系数an=x/x 得到非负的总产出,模型就称为可行的。问为使模 投入系数矩阵A=(ag)nxn 型可行,投入系数应满足什么条件? 产出向量x=(x1,…xn) 需求向量d=(d1,…dn)2 x=(1-4-ld 2)设农业、制造业和服务业的外部需求分别为 3)若三部门的外部需求分别增加1个草位 50,150,100亿元,求三个部门的总产出 求它们的总产出的增量 基本模型x=(1-0-1d记C=(-4 基本模型x=(-4)d 当需求增加△d时,总产出增量△x=C△d 0.150.10.2 A=0.30.050.3 shivan55 0.20.30 MATLAB 5.3. Ink 056341.26760.4930 d=[50150100 若△d=(1,0,0),即农业外部需求增加1单位时,三部门 口x=(139.2801,267.6056208.137) 总产出应分别增加1.3459,0.5634,0.4382单位。 即c的第1列。C的第2,3列给出了什么? (学静学实鉴 (大学数学实验) 4)如果对于任意给定的、非负的外部需求,都能 桃型可行日d20→x≥0午[→0k→可 得到非负的总产出,模型就称为可行的。问为使模 型可行,投入系数应满足什么条件? 矩阵范数定义4口4F 基本模型x=(1-A)d,其中A≥0 4→0=4 横型可行回d≥0→x≥0区-420 4 Ⅱ4l=max∑ (1-AX1+A+…+A)=1-A+ 4<1∑a<1∑x<x(j=1…n) A≥0 A→0(k→叫)日(-4=∑A 若A是由实际敷据得 xy<x,(=,…n)成立 到,只要初始投入非负 问A→0(k→∞)是否成立? 模型可行 8
8 4）如果对于任意给定的、非负的外部需求，都能 得到非负的总产出，模型就称为可行的。问为使模 型可行，投入系数应满足什么条件？ 1）设有n个部门，已知投入系数，给定外部需求， 建立求解各部门总产出的模型。 2）设投入系数如表2所给，如果今年对农业、制造 业和服务业的外部需求分别为50，150，100亿元， 问这三个部门的总产出分别应为多少。 3）如果三个部门的外部需求分别增加1个 单位，它们的总产出应分别增加多少。 实例1 投入产出模型 1）基本模型 投入系数矩阵 ij n n A a = × ( ) 产出向量 T n x ( x , x ) = 1 L 需求向量 T n d ( d , d ) = 1 L ( 1,2, , ) 1 x d xi i n n j ∑ ij + i = = L = 平衡关系 xi : 第i个部门的产出， xij: 第i个部门对第 j个部 门的投入， di : 第i个部门的外部需求 ij ij j 投入系数 a = x / x 产出 投入 部门 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 ( 1,2, , ) 1 a x d xi i n n j ∑ ij j + i = = L = x = Ax + d (I − A) x = d x I A d 1 ( ) − = − x I A d 1 ( ) − 基本模型 = − 2）设农业、制造业和服务业的外部需求分别为 50，150，100亿元，求三个部门的总产出。 ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = 0.2 0.3 0 0.3 0.05 0.3 0.15 0.1 0.2 A T d =[50 150 100] 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个单位， 求它们的总产出的增量。 x I A d 1 ( ) − 基本模型 = − 1 ( ) − 记 C = I − A 当需求增加∆d时，总产出增量 ∆x = C∆d 4）如果对于任意给定的、非负的外部需求，都能 得到非负的总产出，模型就称为可行的。问为使模 型可行，投入系数应满足什么条件？ 模型可行 1 ( )( ) + − + + + = − k k I A I A L A I A ( ) , 0 1 = − ≥ − 基本模型 x I A d 其中A ∀d ≥ 0 ⇒ x ≥ 0 ( ) 0 1 − ≥ − I A A →0 (k →∞) k ∑ ∞ = − − = 0 1 ( ) k k I A A 问 Ak → 0 (k → ∞) 是否成立？ A ≥ 0 1 1 ( 1, ) 1 1 1 A a x x j j n n i ij n i < ⇔ ∑ ij < ⇔ ∑ < = L = = 1 0 0 1 < ⇒ → ⇒ k → k A A A 模型可行 ∀d ≥ 0 ⇒ x ≥ 0 矩阵范数定义 AB ≤ A ⋅ B k k A ≤ A A →0 (k →∞) k ( 1, ) 成立 1 x x j j n n i ∑ ij < = L = 若A是由实际数据得 到，只要初始投入非负 模型可行 ij ij j a = x / x || || max 1 1 ∑= = n i ij j A |a |
实例2一年生植物的繁殖(实验2) 实例2一年生植物的繁殖 ivans 实验2给出的模型为x+Px1+q2 设a1=0.5,a2=0.25,b=0.2,c=10L p=-a,bc, q=-a,b(1-a,)bc 开始有100棵植物,要求50年后有1000棵植物 Ax=B Al=sparse(l: n, I n, p, n, n); A2=sparse( 1: n-1, 2: n, 1, n, n), q p I 0 Fl, n;=[1, 1; s=[-q*x0, -xn] x1=x(1),%输出第2年植物数量101.7097 k=0 n+1; Xx xO, x', xn 求解Ax=B可得第二年(及以后诸年)植物的数量 图4一年生植物繁殖x 实例3汽车刹车距高(续) 布置实验 实验1§2.1建立了刹车距离d与车速v的关系: 数据拟 目的 学会用 MATLAB软件数值求解线性代数方程组,对选 d=k1+k2v2,i=1,2 代法的收做性和解的稳定性作初步分析 2)通过实例学习用线性代数方程组解决简化的实际 车速与刹车距离的实际数据记作(vd,P1,2,,7(见实验1表1) (#)是超定方程组 内率课上布,或参见网络学堂 shivan57 B=,kyy=(41…d 0652
9 实验2给出的模型为 实例2 一年生植物的繁殖(实验2) xk + pxk−1 +qxk−2 =0, k =2,3,L,n 1 21 p =− =− − a bc q a b a bc , (1 ) Ax=B ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = q p q p q p q p p A 1 1 1 1 O O O ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = − − 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 一年生植物的繁殖 1 2 设a a bc = 0.5, 0.25, 0.2, 10 = == 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 xk 图4 一年生植物繁殖xk MATLAB 5.3.lnk shiyan56 实例3 汽车刹车距离(续) 实验1§2.1建立了刹车距离d与车速v的关系： 2 1 2 d = k v + k v d k v k v i n i i i , 1,2, , 2 = 1 + 2 = L 数据拟合 (#) 车速与刹车距离的实际数据记作(vi, di), i=1,2,…,7（见实验1表1) (#)是超定方程组 T n T n n n k k y d d v v v v , ( , ) , ( , ) 1 2 1 2 2 2 1 1 L = = L ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ Φ= × β Φ = β y 1 2 k k =0.6522, 0.0853 = 结果 MATLAB 5.3.lnk shiyan57 布置实验 目的 1)学会用MATLAB软件数值求解线性代数方程组，对迭 代法的收敛性和解的稳定性作初步分析； 2)通过实例学习用线性代数方程组解决简化的实际问题。 内容 课上布置，或参见网络学堂
