4.4龙贝格求积公式.ppt

1、4.4龙贝格求积公式,各节点为,复合梯形(Trapz)公式为,则由上面公式,有,因此有如下递推公式,上式称为递推的梯形公式,递推梯形公式加上一个控制精度,即 可成为自动选取步长的复合梯形公式,思考,外推加速公式,由复合梯形公式的余项公式,可得,复合Simpson公式,因此由复合Simpson公式的余项,可得,即,当然,令,自己 证明,即,当然,同样由复合Cotes公式的余项,得,令,外推 加速 公式,以上整个过程称为Romberg算法,将 上 述 结 果 综 合 后,其中外推加速公式可简化为,Romberg算法的收敛 阶高达m+1的两倍,Romberg算法的代 数精度为m的两倍, Romber

2、g 算法：, ?, ?, ?, ,4.6 高斯求积公式,构造具有2n+1次代数精度的求积公式,将节点 x0 xn 以及系数 A0 An 都作为待定系数。令 f (x) = 1, x, x2, , x2n+1 代入可求解，得到的公式具有2n+1 次代数精度。这样的节点称为Gauss 点，公式称为Gauss 型求积公式。,例：求 的 2 点 Gauss 公式。,代入 f (x) = 1, x, x2, x3,不是线性方程组，不易求解。,证明： “”,x0 xn 为 Gauss 点, 则公式 至少有 2n+1 次代数精度。,对任意次数不大于n 的多项式 Pm(x)， Pm(x) w(x)的次数不大于

3、2n+1，则代入公式应精确成立：,= 0,“” 要证明 x0 xn 为 Gauss 点，即要证公式对任意次数不大于2n+1 的多项式 Pm(x) 精确成立，即证明：,设,求 Gauss 点 求w(x),4 Gaussian Quadrature, 正交多项式族 0, 1, , n, 有性质：任意次数不大于n 的多项式 P(x) 必与n+1 正交。,Step 1：构造正交多项式2,即：,4 Gaussian Quadrature,Step 2：求2 = 0 的 2 个根，即为 Gauss 点 x0 ，x1,Step 3：代入 f (x) = 1, x 以求解 A0 ，A1,解线性方程组，简单。,

4、结果与前一方法相同：, 利用此公式计算 的值, 特殊正交多项式族：,Legendre 多项式族：,满足：,由 有递推,以 Pn+1 的根为节点的求积公式称为高斯-勒让德公式。,例 的两个零点为,则两点高斯-勒让德公式为,计算得到A0=A1=1，因此求积公式为, Gauss 公式的余项：,/* 设P为f 的过x0 xn的插值多项式 */,/*只要P 的阶数不大于2n+1，则下一步等式成立*/,插值多项式的余项,Q：什么样的插值多项式在 x0 xn 上有 2n+1 阶？,A：Hermite 多项式！,满足,Gauss列主元消去法,例1.,用Gauss消去法解线性方程组(用3位十进制浮 点数计算）,

5、解:,本方程组的精度较高的解为,用Gauss消去法求解(用3位十进制浮点数计算),一、Gauss列主元消去法的引入,9999,回代后得到,与精确解相比,该结果相当糟糕,究其原因,在求行乘数时用了很小的数0.0001作除数,主元,如果在求解时将1,2行交换,即,0.9999,回代后得到,这是一个相当不错的结果,例2.,解线性方程组(用8位十进制尾数的浮点数计算),解:,这个方程组和例1一样,若用Gauss消去法计算会有 小数作除数的现象,若采用换行的技巧,则可避免,经过回代后可得,事实上,方程组的准确解为,例:用列主元高斯消去法解线性代数方程组,5.2 Gauss消去法,一、消元与回代计算,对线

6、性方程组,对其增广矩阵施行行初等变换:,定义行乘数,且,定义行乘数,No unique solution exists.,What if we cant find such k ?,No unique solution exists.,定理,若A的所有顺序主子式 /* determinant of leading principal submatrices */ 均不为0，则高斯消元无需换行即可进行到底，得到唯一解。,注：事实上，只要 A 非奇异，即 A1 存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。,二、矩阵的三角分解,1.Gauss消去法消元过程的矩阵描述,行变换

7、相 当于左乘 初等矩阵,由于,令,则,显然若令,则有,因此,从而,故,即,且,顺序主元,定义1. 不带行交换的Gauss 消去法的消元过程,产生 一个单位下三角矩阵L和一个上三角矩阵U,即 该过程称之为,由上述分析不难得到,Hey hasnt GE given me enough headache? Why do I have to know its matrix form?!,When you have to solve the system for different with a fixed A.,Could you be more specific, please?,Factorize

8、 A first, then for every you only have to solve two simple triangular systems and .,Gauss消去法 可以执行,定理1.,在定理中,可能注意到,可能存在,2 Matrix Factorization Matrix Form of G.E.,证明：由1中定理可知，LU 分解存在。下面证明唯一性。,若不唯一，则可设 A = L1U1 = L2U2 ，推出,Upper-triangular,Lower-triangular With diagonal entries 1,注： L 为一般下三角阵而 U 为单位上三角阵

9、的分解称为Crout 分解。 实际上只要考虑 A* 的 LU 分解，即 ，则 即是 A 的 Crout 分解。,例：用矩阵的直接三角分解法（LU分解，L为单位下三角阵、U为上三角阵）解方程组,。,二、Gauss消去法的运算量,计算机作乘除运算所耗时间要远远多于加减运算,且在一个算法中，加减运算和乘除运算次数大体相当,故在衡量一个算法的运算量时只需统计乘除的运算次数,乘法次数：,除法次数：,全部回代过程需作乘除法的总次数为,于是Gauss消去法的乘除法运算总的次数为,数级,Gauss消去法乘除法约为2700次,而如果用Cramer法则的乘除法运算次数约为,或,用行列式定义,用行列式性质,一、直接

10、法概述,直接法是将原方程组化为一个或若干个三角形 方程组的方法，共有若干种,对于线性方程组,其中,系数矩阵,未知量向量,常数项,则,都是三角 形方程组,上述方法称为直接三角形分解法, 道立特分解法 /* Doolittle Factorization */： LU 分解的紧凑格式 /* compact form */,反复计算, 很浪费哦 ,固定 i : 对 j = i, i+1, , n 有,lii = 1,a,固定 j ，对 i = j, j+1, , n 有,b,上述解线性方程组的方法称为 直接三角分解法的 Doolittle法,例1. 用Doolittle法解方程组,解:,由Dooli

11、ttle分解,Doolittle法在计算机上实现是比较容易的,但如果按上述流程运算仍需要较大的存储空间:,因此可按下列方法存储数据:,直接三角分解的Doolittle法可以用以下过程表示:,存储单元(位置),紧凑格式的 Doolittle法,例2. 用紧凑格式的Doolittle法解方程组,解:,所以,定义1.,一、向量和矩阵的范数,显然,并且由于,例.求下列向量的各种常用范数,解:,可以理解为,可以理解为对任何向量范数都成立。, 矩阵范数 /* matrix norms */,(4)* | AB | | A | | B | （相容 /* consistent */ 当 m = n 时),In

12、 general, if we have | AB | | A | | B | , then the 3 norms are said to be consistent.,Oh havent I had enough of new concepts? What do I need the consistency for?,When you have to analyze the error bound of AB imagine you doing it without a consistent matrix norm,例.,不难验证其满足定义2的4个条件,称为Frobenius范数,简称F-

13、范数,而且可以验证,tr为矩阵的迹,类似向量的 2-范数,定义.,简称为从属范数或算子范数,显然，由定义不难推出,算子范数,/* operator norm */,由向量范数 | |p 导出关于矩阵 A Rnn 的 p 范数:,则,矩阵 ATA 的最大 特征根 /* eigenvalue */,例.,求矩阵A的各种常用范数,解:,由于,特征方程为,容易计算,计算较复杂,对矩阵元素的 变化比较敏感,不是从属范数,较少使用,使用最广泛,性质较好, 我们只关心有相容性的范数，算子范数总是相容的。, 即使 A中元素全为实数，其特征根和相应特征向量 /* eigenvector */ 仍可能是复数。将上

14、述定义中绝对值换成复数模均成立。,若不然，则必存在某个向量范数 | |v 使得 对任意A 成立。,Counterexample ?, 谱半径 /* spectral radius */, (A),证明：,由算子范数的相容性，得到,将任意一个特征根 所对应的特征向量 代入,证明：,A对称,若 是 A 的一个特征根，则2 必是 A2 的特征根。,又：对称矩阵的特征根为实数，即 2(A) 为非负实数， 故得证。,对某个 A 的特征根 成立,所以2-范数亦称为谱范数。,证明：, 若不然，则 有非零解，即存在非零向量 使得,线性方程组的误差分析 /* Error Analysis for Linear

15、system of Equations */,求解 时，A 和 的误差对解 有何影响？, 设 A 精确， 有误差 ，得到的解为 ，即,绝对误差放大因子,又,相对误差放大因子,2 Error Analysis for ., 设 精确，A有误差 ，得到的解为 ，即,Wait a minute Who said that ( I + A1 A ) is invertible?,(只要 A充分小，使得,大,显然,即任意方阵的条件数必不小于1,根据算子范数的不同也有不同的条件数:,2 Error Analysis for ., cond (A) 取决于A，与解题方法无关。,常用条件数有：,cond (A

16、)1,cond (A),cond (A)2,特别地，若 A 对称，则,条件数的性质： A可逆，则 cond (A)p 1； A可逆， R 则 cond ( A) = cond (A) ； A正交，则 cond (A)2=1； A可逆，R正交，则 cond (RA)2 = cond (AR)2 = cond (A)2 。,2 Error Analysis for .,精确解为,A1 =,解：考察 A 的特征根,39206 1, 测试病态程度：,此时精确解为,2.0102 200%,cond (H2) =,27,cond (H3) ,748,cond (H6) =,2.9 106,cond (Hn

17、) as n ,注：一般判断矩阵是否病态，并不计算A1，而由经验得出。 行列式很大或很小（如某些行、列近似相关）； 元素间相差大数量级，且无规则； 主元消去过程中出现小主元； 特征值相差大数量级。,三、迭代法的收敛性,设解线性方程组的迭代格式,将两式相减,得,则,因此迭代法收敛的充要条件,可转变为,定理1.迭代格式收敛的充要条件为,4 Convergence of Iterative methods,证明：,“”：对任意非零向量 有,But hey, you dont seriously expect me to compute Bk whenever I want to check the

18、convergence, do you?,根据矩阵与其Jordan标准形及特征值的关系,可知,即,因此,定理.,迭代格式收敛的充要条件为,又因为矩阵的谱半径不超过其任一种算子范数,即,于是又可得到,4 Convergence of Iterative methods,证明：,在计算时,迭代终止的时间可以用上式判别,例.,判别下列方程组用J法和G-S法求解是否收敛,解:,(1) 求Jacobi法的迭代矩阵,因此不能用范数,只能用特征值判断,所以,即Jaobi迭代法收敛,(2) 求Gauss-Seidel法的迭代矩阵,同样用特征值判断,所以Gauss-Seidel迭代法发散,在前面的例子中,G-S

19、法收敛速度比J法要高,但本例却说明G-S法发散时而J法却收敛,因此,不能说G-S法比J法更好,另外,给出系数矩阵对角占优线性方程组的一个结论,定理.,证:,(1)对于Jacobi迭代法,其迭代矩阵为,Jacobi迭代法收敛,(2)对于GS迭代法,其迭代矩阵为,不能使用范数,而用特征值,即,从而,因此,由于,可得,矛盾,GS迭代法收敛,一、简单迭代法(基本迭代法),设线性方程组的一般形式为,依此类推,线性方程组可化为,令,故迭代过程化为,等价线性方程组为,称上式为解线性方程组的Jacobi迭代法(J法),例1.,用Jacobi迭代法求解方程组,误差不超过1e-4,解:,依此类推,得方程组满足精度

20、的解为x12,迭代次数 为12次,x4 = 3.0241 1.9478 0.9205 d = 0.1573 x5 = 3.0003 1.9840 1.0010 d = 0.0914 x6 = 2.9938 2.0000 1.0038 d = 0.0175 x7 = 2.9990 2.0026 1.0031 d = 0.0059 x8 = 3.0002 2.0006 0.9998 d = 0.0040 x9 = 3.0003 1.9999 0.9997 d = 7.3612e-004 x10 = 3.0000 1.9999 0.9999 d = 2.8918e-004 x11 = 3.0000

21、2.0000 1.0000 d = 1.7669e-004 x12 = 3.0000 2.0000 1.0000 d = 3.0647e-005,Jacobi 法和 Gauss - Seidel 法 /* Jacobi 对f (x) 要求不高(只要连续即可) .,无法求复根及偶重根 收敛慢,注：用二分法求根，最好先给出 f (x) 草图以确定根的大概位置。或用搜索程序，将a, b分为若干小区间，对每一个满足 f (ak)f (bk) 0 的区间调用二分法程序，可找出区间a, b内的多个根，且不必要求 f (a)f (b) 0 。,迭代法 /* Fixed-Point Iteration */,

22、f (x) = 0,x = (x),f (x) 的根, (x) 的不动点,思路,从一个初值 x0 出发，计算 x1 = (x0), x2 = (x1), , xk+1 = (xk), 若 收敛，即存在 x* 使得 ，且 连续，则由 可知 x* = (x* )，即x* 是 的不动点，也就是f 的根。,So basically we are done! I cant believe its so simple! Whats the problem?,Oh yeah? Who tells you that the method is convergent?,7.2不动点迭代法(基本迭代法),将非线性

23、方程f(x)=0化为一个同解方程,继续,称上式为求解非线性方程x=(x)的不动点迭代法,则称迭代法收敛,否则称为发散,如果将原方程表示为,与原方程同解,收敛,发散,例.,解:,(1) 将原方程化为等价方程,显然迭代法发散,(2) 如果将原方程化为等价方程,仍取初值,x2 = 0.9644 x3 = 0.9940 x4 = 0.9990 x5 = 0.9998 x6 = 1.0000 x7 = 1.0000,依此类推,得,已经收敛,故原方程的解为,同样的方程 不同的迭代格式 有不同的结果,什么形式的迭代法 能够收敛呢?,迭代函数的构造有关,定理.,(全局收敛性),证:,由条件(1),由根的存在定

24、理,由,由微分中值定理,因0L1，故当k时，序列xk收敛到x*。,证毕.,该定理指出,只要,因此,当,迭代就可以终止,只要构造的迭代函数满足,此时虽收敛但不 一定是唯一根,注：定理条件非必要条件，可将a, b缩小 定义局部收敛性： 设x*为的(x)不动点，若在 x* 的某 邻域 B() = x | | x x* | 有 (x) 连续， 且 | (x*) | 1， 则x0B() 开始的迭代收敛。 即调整初值可得到收敛的结果。,例.,用迭代法求方程的近似解,误差小于10-6,解:,本题迭代函数有两种构造形式,因此采用迭代函数,d1 = 0.1000000 d2 = -0.0105171 d3 =

25、0.115610-2 d4 = -0.126510-3 d5 = 0.139010-4 d6 = -0.150010-5 d7 = 0.100010-6,由于|d7| =0.100010-6 110-6,因此原方程的解为,x7 = 0.090525,x1 = 0.1000000 x2 = 0.0894829 x3 = 0.0906391 x4 = 0.0905126 x5 = 0.0905265 x6 = 0.0905250 x7 = 0.0905251,由定理1的结论可以看出,迭代法收敛就越快,定义.,不可能直接确定,定理.,7.4 牛顿法 /* Newton - Raphson Method */,原理：将非线性方程线性化 Taylor 展开 /* Taylors expansion */,取 x0 x*，将 f (