数值分析[胡良剑]PDF课件.pdf 198页

数值分析 上课 每周三8:15-9:45 二教125 双周三10:05-11:35 图文3号机房 请按机号入座！ 考试成绩 笔试70%，上机考试20%，平时10% 胡良剑 Ljhu@dhu.edu.cn 信息理学楼543 543 数值分析 数值分析是数学与计算机技术结合 的一门学科，研究利用计算机解决 数学问题的理论和方法，是计算数 学的一个重要分支。 zh.wikipedia.org 现代复杂工程技术问题的解决步骤 工程问题 数学模型 设计算法 问题解答 结果分析 上机计算 第一章数值分析的基本概念 §1.1 数值算法的研究对象 §1.2 误差分析的概念 §1.3 数值算法设计的一些要点 例1. 1 (易计算问题) (1) 求解线性方程组AX=B , 其中A 为3阶可逆方 T 阵，X =(x , x , x ) ; 1 2 3 2 * (2) 求代数方程3x +8x−3=0在[0,1]上的根x ； (3) 已知y =P(x)为[x , x ]上的直线，满足P(x )=y , 0 1 0 0 P(x )=y , x ∈(x , x ), 求P(x ) ； 1 1 2 0 1 2 b 1 (4) 计算定积分 I ∫a dx (1ab) ； x (5) 解常微分方程初值问题 y x +y   y (0) 0 解： D (1) Cramer法则 xj j , j 1, 2,3 , 其中D=|A |, D Dj 为由B置换D 的第j 列所得。 * (2)根据求根公式得x =1/3; y 1 −y 0 (3) P(x )= y 0 +x −x (x2 −x0 ) ; 2 1 0 b (4) 根据积分公式得到 I ln ; a (5) 根据线性常微分方程求解公式得 y (x) −x −1+ex 例1.2 (难计算问题) (1) 求解线性方程组AX=B , 其中A 为30阶可逆方 T 阵，X=(x , x ,…, x ) ; 1 2 30 x * (2) 求超越方程xe =1在[0,1]上的根x ； (3) 已知y =f (x)为[x , x ]上的函数，满足f (x )=y , 0 1 0 0 f (x )=y , x ∈(x , x ), 求f (x ) ； 1 1 2 0 1 2 b 1 (4) 计算定积分 I ∫a ln x dx(1ab) ； (5) 解常微分方程初值问题  2 y x =+y   y (0) 0 解：例1.2同例1.1“差不多” ？ (1) 计算量非常大； * (2) 无法求得x 的解析形式，只能求近似值； y 1 −y 0 (3)f (x )≈ y 0 +x −x (x2 −x0 ) 试试； 2 1 0 (4) 无法找到原函数，考虑近似方法； (5) 没有解析解，数值解法求取近似解。 利用计算机！但是……. 计算机的认识能力是有限的 计算机的计算能力也是有限的 可行且高效的算法+计算机！！！ 2 数值算法的特点： 计算机算法 对于给定的问题和设备（计算机），一个算法是 用该设备可理解的语言表示的，对解决这个问题 的一种方法的精确刻画。 计算机算法主要包含数值算法、非数值算法 和软计算方法三类。 三类计算机算法 数值算法主要指与连续数学模型有关的算法， 如数值线性代数、方程求解、数值逼近、数 值微积分、微分方程数值解和最优化计算方 法等；(本课程内容) 非数值算法主要指与离散数学模型有关的算 法，如排序、搜索、分类、图论算法等； 软计算方法是近来发展的不确定性算法的总 称，包括随机模拟、神经网络计算、模糊逻 辑、遗传算法、模拟退火算法和DNA算法等。 数值算法的特点 有穷性 数值性 近似性 数值计算软件 Fortran C++ Matlab §1.2 误差分析的概念 1. 误差限和有效数字 2. 截断误差与收敛性 3. 舍入误差和数值稳定性 4. 数据误差和病态问题 1. 误差限和有效数字 误差和相对误差(定义1.1) * * 设x 是某量的准确值，x 是x 的近似值 * 称δx = x -x 为x 的误差或绝对误差。 * | x -x |≤ε, 称ε为x 的(绝对)误差限或精度， * * δx = (x -x)/x 称为x 的相对误差 r * * |(x - x)/ x |≤ε , 称ε 为x 的相对误差限。 r r 当ε 很小时，ε ≈ε/| x |。 r r 准确位数和有效数字(定义1.2) 设x =±0.a a …a …×10m (m为整数) (1.1) 1 2 n 其中a ~a 为0~9 中一个数字且a ≠0 。 1 n 1 如果 * −k | x -x |≤0.5×10 (1.2) 即x 的误差不超过10-k位的半个单位 则称近似数x准确到第k位小数，并说x有m+k位 有效数字。 例1.3 (误差限和有效数字) 圆周率π=3.1415926…。 x =3.14; x =3.141; x =3.142; x =3.1414 1 2 3 4 解(1) π−x =0.15926…×10−2, |π−x | ≤ 0.5×10−2, 1 1 有3位有效数字； (2) π−x =0.5926…×10−3, |π−x |≤0.5×10−2,有3位 2 2 有效数字； (3) π−x = −0. 4073…×10−3, |π−x |≤0.5×10−3,有4 3 3 位有效数字； (4) π−x =0.1926…×10−3, |π−x |≤0.5×10−3, 有4位 4 1 有效数字。 有效数字概念的通俗定义 * * 设x 是某量的准确值，x 是x 的近似值，如果 在从第一个非零数字开始的第n位进行四舍 * 五入(即考虑第n+1位是舍还是入？)，x 和x 的结果完全一致，则称x有n位有效数字。 与定义1.2的区别 2. 截断误差与收敛性 截断误差：一个无限的数学极限过程用有限次运算 近似计算产生的误差。 例（无限） 2 n x x x e 1=+x + ++ + 2! n! 近似计算（有限） 2 n x x x ≈ ( ) 1+ + ++ e S x x n 2! n! 截断误差（余项公式） n+1 x ex S x R x eξ − ( ) ( ) n n ξ在0与x之间 (n +1)! 算法的收敛性：该算法总可以通过提高计算 量(增加项数)使得截断误差任意小。即 余项→0 3. 舍入误差和数值稳定性 舍入误差: 由于机器字长的限制而产生的误差 机器数(二进制0-1,离散) 规格化浮点式: 阶码m,字长t ,尾数β (β =1) 1 m ±2 ×0.β β …β , m=±α α …α 1 2 t 1 2 s 0 单精度32位(4字节): t=23,s=7, 符号2位, 表示范围± −39 38 -128 128 -23 2.9×10 ∼3.4×10 (即2 ∼ 2 ),识别差2 =1.19 ×10−7 双精度64位(8字节) : t=52,s=10, 符号2位,表示范围 ± 5.56×10−309∼1.79×10308 (即2-1024 ∼ 21024),识别差 -52 −16 2 =2.22 ×10 上溢出±∞和下溢出0 数值稳定性 误差传播问题: 设函数y =f (x , x , …, x )是一 1 2 n *  x x 个算法或模型， 是变量x 的准确值，而 i i i 是变量x 的近似值。如果 |  | | *  | , δx x =−x ≤ε i i i i i 且f 的计算过程中没有新的误差产生，那么计     算结果 y f (x ,x , ,x ) 具有怎样的精度？即 1 2 n  * * *    |δy | | f (x ,x ,,xn ) =−f (x ,x ,,xn ) |≤? 1 2 1 2 算法的数值稳定: 计算过程中舍入误差 不会被严重放大 误差的传播 线性情形用严格估计 n  * * *    |δ y | | f (x ,x ,,x ) =−f (x ,x ,,x ) |≤ |a | ε 1 2 n 1 2 n ∑ i i i 1 非线性情形用线性近似 * n f n ∂ f x * ∂ *    i  δ y ≈∑ ( ) δxi δr (y ) ≈∑( ) * δr (xi ) i 1 ∂xi i 1 ∂xi y ∂f * 绝对误差传播主要取决于条件数 | ( ) | ∂x i * ∂f * x 相对传播主要取决于条件数 | ( ) i | * ∂xi y 四则运算 δ(a±b)=δa±δb, δ(a±b)= [a/(a±b)]δa+[b/(a±b)]δb(近似数相减 r r r 不稳定) δ(ab) ≈bδa+aδb δ(ab) ≈δa+δb r r r 2 δ(a/b) ≈(1/b)δa−(a/b )δb (b ≈0不稳定) δ(a/b) ≈δa−δb r r r 例1.5 (数值稳定性) 1 I n ∫x n ex −1dx n=0, 1, …, 20 0 1 1 估计 ≤I n ≤ ( 1) 1 n + e n + 算法一：递推公式 I =1−nI , n=1,…,20 n n-1 I =1-1/e →I →I →… →I 误差很大 0 1 2 20 δ = −nδ , δ =(20!) δ ，不稳定 n n-1 20 0 算法二：递推公式 I =(1−I )/n , n= 20 ,…,1 n-1 n I20≈估计式中点→I19 →… →I1 →I0 误差很小 δ = −δ /n, δ =δ /(20!) ，稳定 n-1 n 0 20 4. 数据误差和病态问题 例1.6 (病态问题)  49 3 21 949 x + x + x  36 1 4 2 40 3 360 1.3611x +0.7500x +0.5250x 2.6361  3 36 3 1061  1 2 3 x + x + x 0.7500x +0.4360x +0.3000x 1.4736  1 2 3 1 2 3  4 144 10 720  0.5250x +0.3000x +0.2136x 1.0386 21 3 769 3739  1 2 3 x + x + x  1 2 3 40 10 3600 3600 x =x =x =1 x =1.1788, x =-0.0006, x = −1.9658. 1 2 3 1 2 3 病态问题: 很小的变化数据却导致解产生了很大的 变化。 区别：收敛性和数值稳定性主要源于算法，病态性 主要是模型本身的原因 。 §1.3 数值算法设计的一些要点 设计算法基本原则 精度：收敛性、稳定性 速度：计算量、存储量，收敛速度 注意事项 病态问题 速度细节(加法、乘法，函数) 稳定性细节(近似数相减，分母接近0，大数吃小 数) 例1.8(相近的数相减) x=π×10-8 ， 1 1−x → 2x 2 y − y 1+2x 1+x (1+2x)(1+x) 例1.9 (大数吃小数) 1000 y =123456+ ∑xi , 0.01≤x ≤0.04 i i 1 存储细节(降维) 死循环 实数相等比较 2 2 2 2 sin x=1-cos x →|sin x-(1-cos x)|≤ε 中间结果 作业 P11, ex1, ex2, ex5, ex7 第二章 数值代数 §2.1 Gauss消去法 §2.2 直接三角分解法 §2.3 范数和误差分析 引例：Cramer法则不可行 a x ++a x b  11 1 1n n 1    a x ++a x b n1 1 nn n n Cramer法则 D j , 1, , xj j n D n20时，计算量太大，现实上不可行 Cramer法则数学上很重要，计算上无价值 §2.1 Gauss消去法 1 理论基础 2 顺序Gauss消去法 3 选主元技术 4 追赶法 1. 理论基础 引理2.1 a x ++a x b  11 1 1n n 1    a x ++a x b m1 1 mn n m 行初等变换   解 同 (A,b ) →(A,b ) ~ ~  ~ a x + +a x b  11 1 1n n 1   ~ ~ ~  a x  a x b + +  m1 1 mn n m 证明(P14) 2 顺序Gauss消去法 例2.1 x1 + x 2 +x3 6  −x +3x +x 4 1 2 3  2x −6x +x =−5  1 2 3 消元过程 回代过程 顺序Gauss消去法 b a(1) x ++a(1) x a(1)  11 1 1n n 1,n+1    (1) (1) (1) a x ++a x a n1 1 nn n n ,n+1 Ax 消元过程(第1步) 用矩阵初等行变换化系数矩阵为三角形 a(1) a(1) a(1)  a(1) a(1)   11 12 13 1n 1,n+1  (1) (1) (1) (1) (1) a21 a22 a23  a2n a2,n+1  a(1) a(1) a(1)  a(1) a(1)   31 32 33 3n 3,n+1          (1) (1) (1) (1) (1)  a a a  a a , +1  n1 n 2 n 3 nn n n  (1) a (2) (1) (1) i1 , , a a =−l a l ij ij i1 1 j i1 (1) a 11 i 2, , n; j 2, , n =+1 消元过程(第2步) a(1) a(1) a(1)  a(1) a(1)   11 12 13 1n 1,n+1  (2) (2) (2) (2)  0 a22 a23  a2n a2,n+1  → 0 a(2) a(2)  a(2) a(2)   32 33 3n 3,n+1          (2) (2) (2) (2)  0 a a  a a n n+  n 2 n3 nn , 1  (2) a (3) (2) (2) i2 , , a a =−l a l ij ij i2 2 j i2 (2) a 22 i 3, , n; j 3, , n =+1 消元过程(第k步)编程用计算公式 (k ) a (k +1) a (k ) −l a (k ) l aik ij ij ik kj ik (k ) a kk i k +1, , n;j k +1, , n +1; k 1, , n −1 n−1 乘法次数 − − + ∑(n k )(n k 1) k 1 n−1 除法次数∑(n −k ) k 1 消元过程(结果) a(1) a(1) a(1)  a(1) a(1)   11 12 13 1n 1,n+1  (2) (2) (2) (2)  0 a22 a23  a2n a2,n+1  上三角矩阵 → 0 0 a(3)  a(3) a(3)   33 3n 3,n+1          (n) (n )  0 0 0  a a  nn n ,n+1  主元素a(k ) ≠0，k 1,2, , n kk 上三角方程组  (1) (1)  (1) (1) a x +a x + +a x a 11 1 12 2 1n n 1,n+1 (2) (2) (2)   a22 x2 + +a2n xn a2,n+1     a(n )x a (n ) 回代 nn n n n+  , 1 回代过程  (n) (n) x a a , 1  n n n+ nn n (k ) (k ) (k ) ⇒xk (ak n − ∑akj xj ) a , +1 kk  j k +1 (k n −1, ,1)  n−1 乘法次数 (n −k ) ∑ k 1 n−1 除法次数 1 ∑ k 1 计算量(乘除法次数) 消元 n−1 n(n −1)(2n +5) ∑(n −k )(n −k +2) k 1 6 回代 n n n n + ( 1) n k i ∑( − +1) ∑ k 1 i 1 2 总和 计算量主要在消元 n3 2 n n3 比较Cramer法则 +n − ≈ 3 3 3 Gauss 法快很多 可行性 如果A=(a ) 的顺序主子式均不为零，顺序 ij n×n Gauss消去法求解可行。 3 选主元技术 为什么要选主元素？ 主元素=0，计算就不能进行。 主元素≈0, 计算过程数值不稳定。 怎样选主元素？ 把绝对值大的数调到对角线上 例2.2 3位有效数字 0.001x +x 1  1 2  x +x 2  1 2 例2.3(选列主元素Gauss消去法)  1 1 1 6  r ↔r  2 −6 1 −5   1 3 −1 3 1 4    −−−−−−→ −1 3 1 4      2 −6 1 −5 (选主元)  1 1 1 6  r +0.5r 2 1 2 −6 1 −5  2 −6 1 −5  r −0.5r r ↔r 3 1   2 3   −−−→ 0 0 1.5 1.5 −−−→ 0 4 0.5 8.5     0 4 0.5 8.5  0 0 1.5 1.5  可行性 如果A=(a ) 的行列式不为零，选列主元素 ij n×n Gauss消去法求解可行。 比较： 如果A=(a ) 的顺序主子式均不为零， ij n×n 顺序Gauss消去法求解可行。 4 追赶法 三对角方程组 b x +c x d 1 1 1 2 1  a x +b x +c x d 2 1 2 2 2 3 2       a x +b x +c x d − − − − − −  n 1 n 2 n 1 n 1 n 1 n n 1  a x − +b x d  n n 1 n n n 顺序Gauss消去法应用于三对角线性方程组 得到所谓追赶法，其中消元过程为 “追”， 回代过程为“赶”。 4 追赶法 ~ 追 b b , ~ 1 1 d 1 d 1 a ~ ~ ~ l k b b −l c k ~ , k k k k −1, d d −l d b k k k k −1 k −1 k=2, …, n 赶 ~ ~ d d −c x n k k k +1 x ~ x n k ~ b b n k k= n-1, …, 1 追赶法不对零元素计算，只有2(n-1)+(n-1)+n+(n-1) =5n-4次乘除法计算量 注意:追赶法假定主元不为0, 计算中不选主元. §2.2 直接三角分解法 高斯消去法的矩阵表示 LU分解法 平方根法 改进的平方根法 高斯消去法的矩阵表示 对一个矩阵施行一次行变换，相当于左乘一个相 应的初等矩阵。 顺序高斯消去法相当于矩阵三角分解A=LU P21 例题 P P P A=U 3 2 1 选列主元高斯消去法相当于三角分解PA=LU, 其 中P为行置换矩阵(详见课本38页) LU分解 A=LU (Doolittle分解)  1  u11 u12  u1n      L l21 1  U  u22  u2n           l l  1 u  n1 n 2   nn  解方程组 Ly b Ax b =⇒LUx b =⇒ Ux y  P22 例2.5待定系数法 LU分解计算顺序 待定系数法(依次显式计算，无须解方程组) a11 a12  a1n   1 u11 u12  u1n       a a a l 1 u  u  12 22 1n  21  22 2n                   a a  a l l  1 u  n1 n 2 nn   n1 n 2  nn  LU分解法vs. Gauss消去法 LU分解法和Gauss消去法具有相同的可行性条 件, 基本相同的计算量和计算精度。 LU分解法具有比Gauss消去法更好的设计灵活 性。 当多次求解具有相同系数矩阵和不同右段向量的线 性方程组 ； 防止“大数吃小数”； LU分解存储可使用紧凑格式 。 三角分解的其他形式: Crout分解等. 平方根法 l11    l l A=LL T(Cholesky分解) L  21 22         l l  l 基本方法: 待定系数法  n1 n 2 nn  手算：由顺序主子式先求对角线 D lkk k ,k 1,, n D 编程: 使用公式 k −1 k −1 k −1 1 l a − l 2 lik (a −∑l l ) kk kk ∑kj l ik ij kj j 1 kk j 1 i=k+1, …n, k=1,…,n 乘除法次数n(n+4)(n-1)/6, 开方n次 P24 例2.6 数值稳定性 可行性条件: 对称正定. 直接验证知：在平方根法中， k a l 2 ⇒l ≤ a ，k=1, …,n, j =1, …, k kk ∑kj | kj | kk j 1 中间量lkj不会出现放大, 从而平方根法是数值 稳定的。不必选主元. 改进的平方根法 T  1  A=LDL ,   d1  L l21 1  D        ,  dn  l l  1  n1 n 2  待定系数法