数值计算方法与算法

1、数值计算方法与算法教学目标 掌握常用的数值计算方法 掌握计算方法的数学原理 学会选择恰当的计算方法 学会使用流行的计算软件教学计划时间：7:50-8:35地点：21062-27 绪论3-06多项式插值3-13 多项式插值3-20分段和样条插值3-27 数值微分4-03数值积分4-10 数值积分4-17最小二乘法4-24 方程求根5-01五一放假5-08 方程求根5-15求解线性方程组5-22 求解线性方程组5-29计算矩阵特征值6-05 计算矩阵特征值6-12微分方程数值解6-19 微分方程数值解6-26复习7-04 期末考试7-13报送成绩教材及参考 数值计算方法与算法，张韵华、奚梅成、陈效

2、群编，科学出版社，2006。 科学计算导论，Michael T. Heath著，张威、贺华、冷爱萍译，清华大学出版社， 2005。 数值计算方法解题指导，张韵华编，科学出版社，2003。 网络教学资源 http:/联系方式教师：王新茂 理化大楼16#016，3607565助教：杨荣136-15607693第0章 绪论 数学建模 数值计算实际问题数学问题近似解 什么是数值计算方法？ 什么是“好的”数值计算方法？误差小 误差分析耗时少 复杂度分析抗干扰 稳定性分析 误差的类型绝对误差真实值近似值相对误差绝对误差真实值 误差的来源原始误差、截断误差、舍入误差输入计算输出真实值近似值xxxxf)(xf

3、y fffyxfy)( 一些例子：计算地球的体积计算计算 如何减小计算误差？选择好的算法、提高计算精度 范数的定义满足非负性,齐次性,三角不等式的实函数334RV 715131143223333)(),(yxyyxxyxyxf 常用的向量范数 常用的矩阵范数 矩阵的谱半径 例：计算矩阵 的范数和谱半径。 例：范数在误差估计中的应用pxxxppnpp111，pxAxAppp1sup，nA,max)(14321A作业一、145页习题6第1,2题.作业二、利用公式 编写两个计算ex的C程序(一个用单精度, 另一个用双精度).令x=1,5,10,15,20, 比较它们和库函数exp(x)之间的运行时间

4、和计算误差. 思考如何减小误差？0!nnxnxe第1章 插值函数逼近用未知函数f(x)的值构造近似函数(x)。要求误差小、形式简单、容易计算。常用的函数逼近方法 插值：(xi)=yi, i=0,1,n. 拟合：|(x)-f(x)|尽可能小通常取 (x) = a00(x) + + ann(x)，其中i(x)为一组基函数。多项式插值 给定平面上n+1个插值点(xi,yi), 构造n次多项式(x), 满足(xi)=yi, i=0,1,n.nhxnhxnnaaaxxaxaax)()()(1010，或01110101100)()(111axaxxaxaxyyyaaaxxxxxxnnnnnnnnn单项式插

5、值)()()()(),()()()(101100ininnxxxxxxxxxLxLbxLbxLbxLagrange插值nnnniiiiiiiinnnnxxbxxbxxxxxxxxxxxxxybyyybbbxLxLxL00011010101100)()()()()()()()()(001111010111)()()()()()(1)(11cxxxxcxxcxyyycccxNxNxNnnnnnnnnn)()()(),()()(110110iinnxxxxxxxNxNcxNccxNewton插值差商表012n012.n0 x1x2xnx0101xxyy1212xxyy11nnnnxxyy01, 1,

6、 1xxnnnnn021 , 12, 1xx 21, 1, 1nnnnxx,jijijxxfk阶差商0110110,xxxxxfxxxfxxfkkkkkny2y1y0y差商的性质 以x0,xn为节点的n次插值多项式(x)的首项系数等于fx0,xn。证明：分别以x0,xn-1和x1,xn为节点构造n-1次插值多项式1(x)和 2(x)，则有对n用归纳法。 fx0,xn与x0,xn的顺序无关。)()()(20010 xxxxxxxxxxxnnn误差估计：证明：设，则有n+2个零点。根据中值定理，存在于是。)()()!1()()()()(0)1(nnxxxxnfxxfxR)()()(0nxxxxxK

7、xR)()()()(0nxtxtxKtftgbagn，0)()1()!1()()()1(nfxKnHermite插值 给定平面上n+1个插值点(xi,yi,mi), 构造2n+1次多项式(x), 满足(xi)=yi, (xi)=mi, i=0,1,n.nnnnnnnnnnnnnnmmyyaaaaxnxxnxxxxxxxxaxaax0012210220012212020012110) 12(210) 12(21011)(单项式基函数Lagrange基函数)(,)(2)()()()(2302iiiiiiiiiiiniiiiixLycxLyxLmbxLcxxbx误差估计：证明：设 ，则有2n+3个零

8、点。根据中值定理，存在于是。220)22()()()!22()()()()(nnxxxxnfxxfxR220)()()(nxxxxxKxR220)()()()(nxtxtxKtftgbagn，0)()22()!22()()()22(nfxKnRunge现象：并非插值点取得越多越好。解决办法：分段插值-1-0.50.510.511.52三次样条插值 给定平面上n+1个插值点(xi,yi), 构造分段三次多项式(x), 满足(xi)=yi, (x)可微，”(x)连续。作业一、习题1第2,4,6,8,10,12,14,16题。作业二、在半圆 上随机选取10个点, 构造插值多项式, 画出函数图像, 并

9、比较3种插值方法的计算误差。作业三、思考3种插值系数之间的关系。比较3种插值方法的优缺点和应用范围。21xy第2章 数值微分和数值积分数值微分 差商法向前差商向后差商中心差商1212)()()(xxxfxfxfhxfhxf)()(hhxfxf)()(hhxfhxf2)()( 插值法在 x 附近取点(xi,f(xi)构造插值多项式。 样条法在 x 附近取点(xi,f(xi)构造样条函数。 f(x)(x)niijjijjijiniijjijixxxxxxxfxxxxxxfx001)()()()(例：用中心差商公式计算f(xi)。例：用向后差商公式计算f(0.2), f(0.4)。x0.00.10.

10、20.30.4f(x)2.01.9f(x)f”(x)x0.00.4f(x)0.8187310.9048371.0000001.1051711.221403f(x)例：设xi=x0+i*h, i=1,.,n。计算(xk)。解：kiikikjkkiijjikjijjkikjjkkkniijjijjijiiniknkikhxfjkhxfxxxxxfxxxfxxxxxxxxfx)!( !)!( !) 1()(1)()()()(1)()(1)()(,0 误差估计前后差商中心差商插值微分26)()(2)()(hfxfhhxfhxf hfxfhxfhxf2)()()()(

11、 ijjiniinxxnfxfxxxxxxKxfx)()!1()()()()()()()()1(0数值积分 插值法 niibaijjijbaniiijjijxfdxxxxxdxxxfxxxxx00)()()()()()()(00nnbaxfxfdxxf 若积分公式对任意m次多项式都取等号，则称积分公式具有至少m阶的代数精度。 插值型积分公式的代数精度n。 当积分节点 x0,.,xn 给定时，代数精度n的积分公式唯一。mnmnmncccxxxx1000011例：设xi=a+i*h, i=0,.,n, h=(b-a)/n。计算Newton-Cotes积分解：nhnnnnnhnhnhnhhnhhnn

12、nnnnnnnn11211101112211011010111)()()(00111badxx)(特别，当n=1,2时，积分公式分别称为梯形公式Simpson公式2)()()(1bfafabI6)()(4)()(22bffafabIbana1a2a3a4a5121/64/61/631/83/83/81/847/9032/9012/9032/907/90 误差估计特别，梯形公式和Simpson公式的误差为代数精度1代数精度3even is )()()!2()(odd is )()()!1()()()()()(0)2(0)1(0nxdxxxxxnfndxxxxxnfdxxxxxxKdxxdxxfb

13、annbannbanbaba5)4(231)(2880)()(12)(abfEabfE 复化数值积分mIIbadxxfdxxfdxxf)()()(1 梯形公式 Simpson公式32101)(12)(2)()()(abnfxfxfnabdxxfniiiba 54)4(1022122)(2880)(6)()(4)(2)(abnfxfxfxfnabdxxfniiiiba Richardson外推法我们要计算假设则有比 和 更高的精度。 误差估计)(lim0hFah)()(qphOhbahF)(1)()(qpphOahFhF)(hF)()()(1)(qpphOhFhFahF)( hF Romberg

14、积分公式 等分的梯形公式， 瑕积分 重积分 Gauss-Legendre积分kkR21 , 3 , 21411, 11,jRRRjjkjkjk，, 2 , 1kdxxf)( dcbadydxyxf),(定理：假设 满足则插值积分公式具有2n+1阶的代数精度。 证明：课本21页性质1.3：若f(x)为m次多项式，则fx0,.,xn,x为m-n-1次多项式。bannbabadxxxxxxxxfdxxdxxf)()(,)()(00)()()(0nxxxxxgnidxxgxbai, 1 , 00)(， 求多项式空间在内积下的标准正交基。解法1：对任意基作Gram-Schmidt正交化。解法2：对任意度

15、量方阵作相合对角化。解法3：求解正交关系的线性方程组。解法4：Legendre多项式badxxgxfgf)()(),()()()( !1)(nnnnnnbxaxdxdabnxL作业一、习题2第111题。作业二、使用各种方法计算的值，其中 0 x1，要求误差1e-9。01)(dtetxtx第3章 曲线拟合的最小二乘法曲线拟合 对区间 I 上的连续函数 f ，构造特定类型的函数 使 f 。 对离散数据序列 (xi, yi), i=1,2,m，构造特定类型的函数 使 (xi)yi 。最小二乘法 求 使最小。 求 使最小。Idxxfx2)()(2)(iiiyx多项式拟合其中 是标准正交基，。 求使 最

16、小。)()()()(1100 xaxaxaxnnIiidxxxfa)()()(xinnxaxaax10)(),(10naaa2YVmnmmnnyyyYxxxxxxV212211111，YVPOQVQOPVrr111111奇异值分解Moore-Penrose广义逆矛盾方程组的解YV其他类型的离散数据拟合 xbax)()cos(cos)(10nxaxaaxnnnmmxbxbbxaxaax1010)(作业一、习题3第15题。作业二、对下列数据，用最小二乘法求形如的经验公式。 dcbxax)cos()(x0.40.5y0.90590.56320.38350.42440.6730 x0

17、.1.0y1.04901.43151.69731.76081.6016第4章 非线性方程求根问题 求f(x)=0在区间a,b内的实根 求f(x)=0在x0附近的一个实根 求f(x)=0在x0附近的一个复根 求多项式f(x)=0的所有复根 求非线性方程组的根方法 用近似函数(x)的根逼近f(x)的根。 二分法已知f(a)f(b)0，设c=(a+b)/2。若f(a)f(c)0则根在b,c内。当|f(c)|或|b-a|时，输出c。迭代步数：O(log2) 不动点当| (x)|L1时，|xk+1-|L|xk-|。当|xn+1-xn|时，输出 xn。迭代步数：O(logL)(1kkx

18、xLipschitz常数线性收敛 Newton法（一阶Taylor展开）211)()(2)()()()()( kkkkkkkxxffxxfxfxxxfxf当|f(xk)|或|xk+1-xk|时，输出xk+1。迭代步数：O(loglog)二次收敛 Newton法（p重根情形）)()()()()()()()()()()()()()()(211kkkkkkkkkkkpxxgxgpxxgxgxxfxfpxxxgxgxpxfxfxgxxf用Newton迭代法求 f(z)=z32z+2 的根。当初值分别位于红、蓝、绿色区域时，迭代收敛到三个根。当初值位于黑色区域时，迭代陷入死循环010。图片引自John

19、Hubbard, Dierk Schleicher, Scott Sutherland, How to find all roots of complex polynomials, Inventiones mathematicae 146, 1-33 (2001). 弦截法（线性插值）)()()(111kkkkkkkxfxfxfxxxx当|f(xk)|或|xk+1-xk|时，输出xk+1。迭代步数：O(loglog) 弦截法的收敛速度618. 12151)|(|)(2)()(lim)()(2)()(2)()(211111111111 rrrxOxffxxxxxxxffxxxxxfxfriikk

20、kkkkkkkkkkkk，则有假设 Newton法解非线性方程组 求的所有复根等价于求 x1,xn 使 f(t)=(t-x1)(t-xn)。OXF)()(|)(12XFXXOXXFjijiXFXFnnnttaatf110)()()(jiijiiixxxfxx 其他求根方法Brent（反插值 x=(y)）Halley（二阶Taylor展开）Muller（二次插值）有理插值作业一、习题4第1、3、5、7、8题。作业二、比较各种求根方法的优缺点。第5章 解线性方程组的直接法问题：求解n元线性方程组AX=B。方法？速度？精度 ？存储？ 下三角方程组 上三角方程组n(n-1)/2次加减法，n(n+1)/

21、2次乘除法。. 1 ,/ )(,11,niaxaxabxiinniiiiii，., 1/ )(11,11niaxaxabxiiiiiiii， Gauss消元法解一般方程组(2n3+3n2-5n)/6次加减法，(n3+3n2-n)/3次乘除法。nnnnnnnnnnnnnnnnnnnnnbababababaabaaabaabaabaaabaaabaaabaaa0022211122221112112222211121121222221111211 追赶法解三对角方程组3n-3次加减法，5n-4次乘除法。nnnnnnnnnnnnndadadadabdadbadacbdadbadacbdacdba022

22、111221111221111222111线性方程组解的精度矩阵条件数bbAAOAAbAbAbbAAbAbAAAObAbbAAAAAOAAAAAAIAAAIAA1111111111111111111)()()()()()()()(1)cond(AAAGauss消元法的实质是LU分解 存在性？A的顺序主子式0。 唯一性？L1U1=L2U2 L1-1L2=U1U2-1 对角 精确度？A-1b的相对误差(L,U,b)的相对误差cond(L)cond(U)。nnnnnnnnuuuuuullllllULA2221121121222111nnnnnnuuuuuulllA222112112121111111

23、211221222111nnnnnnuuullllllADolittle分解Courant分解1, 1 ,1111112112212121ijijnnnnnulQPuuudddlllPAQ是置换阵，全/列/行主元分解LDLT分解、Cholesky分解1111112112212121nnnnnllldddlllAQR分解jjijnnnnnnnnnnnnnnnnnnnnrrQPrrrrrrqqqqqqqqqAPQrrrrrrqqqqqqqqqA是正交阵是置换阵是正交阵2221121121222211121122211211212222111211SVD分解Givens旋转Householder反射

24、1cossinsincos1QTIQ22是正交阵VUVnUA,1作业一、习题5第1-8题(1)。作业二、计算100阶Hilbert矩阵H的逆矩阵A以及AH-I的元素平方和。199110111001101131211001211H第6章 解线性方程组的迭代法 Jacobi迭代 Gauss-Seidel迭代njjnjnnnnjjjjjjnnnnnnnnnnxabxaxabxaxabxabxaxaxabxaxaxabxaxaxa22222211111122112222212111212111ijjijiiiixabax1ijjijijjijiiiixaxabax1 迭代法解线性方程组 A X = B

25、A Xk+1 B = C (A Xk B)C 称为 Conditioner，满足 (C)1或|C|1通常取 C = I A -1，其中 A，于是Xk+1 = Xk -1 (A Xk B) Jacobi迭代： = D定理：A行对角优、或A列对角优 Jacobi迭代收敛。 Gauss-Seidel迭代： = D + L定理：A行对角优、或A列对角优、或A正定 Gauss-Seidel迭代收敛。 松弛迭代： = w-1D + L定理：松弛迭代收敛 0w2定理：A正定且0w2 松弛迭代收敛 Newton迭代求 A-1：Xk+1 = 2 Xk Xk A Xk作业一、145页习题3、6、7。作业二、用迭代

26、法计算10阶Hilbert矩阵H的逆矩阵A以及AH-I的元素平方和。第7章 计算矩阵的特征值和特征向量问题1：求复方阵的模最大(最小)特征值。方法：幂法、反幂法问题2：求复方阵的所有特征值。方法：QR 迭代问题3：求Hermite方阵的所有特征值。方法：Jacobi 方法幂法 当 A 只有一个模最大的特征值max ，并且 x0 与max 的特征向量 amax 不正交时 当 A 的模最大的特征值都相同时，以上迭代仍然收敛。max1maxkkkkkkHkkHkkrAxrAxxxxAxxr， 当 A 的模最大的特征值各不相同时，可以选取数 s 使 A - s I 的模最大的特征值只有一个。 当 A

27、恰有 m 个模最大的特征值时，有 R 的特征值就是 A 的模最大的特征值。XRAXXRAXXXXRkkkkHkkHkk111)()(反幂法 当 A 只有一个模最小的特征值min ，并且 x0 与min 的特征向量 amin 不正交时 计算 A - s I 的模最小的特征值等价于计算 A 的最接近 s 的特征值。min111min1kkkkkkHkkHkkrxArxAxxxxAxr，QR 迭代 利用 QR 分解，酉相似 A 为上三角。 QR 迭代的本质是幂法 当 时，QR 迭代收敛。 可对 A - s I 作 QR 分解，加速收敛。kkkkkkkkkQRQAQARQA1111RRQQAkkkn2

28、1Jacobi 方法 通过 Givens 旋转，逐渐减小非对角元。本质是 2 阶 Hermite 方阵的酉相似。 Jacobi 方法具有 2 阶收敛速度。*00*00*cossinsincoscossinsincosuvvuaaaauvvuaaaajjjiijiijjjiijii复矩阵的奇异值分解 A = UV 一般方法AH A = VH2 V 或 A AH = U2 UH QR 迭代 Jacobi 方法计算 2 阶方阵的 SVDTkkkkkkkTkkkkRQAQARQRRQA,11作业一、167页习题1(3)、2(2)、3(4)。作业二、计算10阶Hilbert矩阵的正交相似标准形以及过渡矩

29、阵。第8章 常微分方程数值解问题：求解一阶常微分方程的初值问题：解法：化微分方程为积分方程bxayayyxfxy0)(),()(1)(,()()(1nnxxnndttytfxyxyEuler折线法 向前Euler公式 向后Euler公式Picard迭代 中心Euler公式 梯形公式改进的Euler公式),(1nnnnyxfhyy),(111nnnnyxfhyy),(211nnnnyxfhyy),(),(11211nnnnnnyxfyxfhyy),()1(11)(1knnnknyxfhyy),(),(),()0(1121)1(1)0(1nnnnnnnnnnyxfyxfhyyyxfhyyRunge