计算方法教案

第1章误差分析与数值计算3§1.1引言3§1.2绝对误差与相对误差、有效数字9§1.3近似数的简单算术运算12§1.4数值计算中误差分析的一些原则13第2章非线性方程（组）的近似解法15§2.1引言15§2.2根的隔离16§2.3对分法16§2.3对分法17§2.4迭代法19§2.6弦截法21§2.6弦截法22§1.7用牛顿法解方程组23本章小结25第3章线性方程组的解法26§3.1引言26§3.2高斯消去法28§3.3矩阵的LU分解31§3.4对称矩阵的LDLT分解32§3.5线性方程组解的可靠性33§3.6简单迭代法34本章小结43第4章矩阵特征值与特征向量的计算44§4.1引言44§4.2幂法和反幂法45§4.3雅可比方法46§4.4QR方法*51本章小结52第5章插值与拟合53§5.1引言53§5.2插值多项式的存在和唯一性54§5.3拉格朗日插值多项式55§5.4均差插值公式57§5.5差分等距结点插值公式59§5.6爱尔米特插值公式61§5.7分段低次插值62§5.8三次样条函数63§5.9曲线拟合的最小二乘法67本章小结70第6章数值积分和数值微分71§6.1引言71§6.2牛顿一科特斯型积分公式72§6.3复合求积公式74§6.4龙贝格求积公式77§6.5高斯求积公式78§6.6二重积分的数值积分法80§6.7数值微分81本章小结83第7章常微分方程的数值解法84§7.1引言84§7.2欧拉法和改进的欧拉法85§7.3龙格-库塔方法86§7.4线性多步法89§7.5算法的稳定性与收敛性91§7.6微分方程组和高阶微分方程解法92本章小结94第1章误差分析与数值计算§1.1引言1、课程任务和目的：在第七届国际软件工程学术会议上，“计算方法”被列入应用方法学的研究领域，强调了计算方法的研究应用与软件方法学的研究密切结合。这就说明了计算方法与软件之间的联系以及在应用软件研制中的地位与作用，计算方法是研究各种数学问题求解的数值计算方法。在计算机成为数值计算的主要工具的今天，则要求研究适合于计算机使用的数值计算方法。计算方法就是研究用计算机解决数学问题的数值方法及其理论，它的内容包括函数的数值逼近、数值微分与数值积分、非线性方程值解、线性方程组数值解、常微和偏微数值解等，即都是以数学问题为研究对象的。因此，计算方法是数学的一个分支，只是它不象纯数学那样只研究数学本身的理论，是把理论与计算紧密结合，着重研究数学问题的数值方法及其理论，计算方法是计算机应用和软件研制开发的重要组成部分，通过本课程的学习和上机实习，使学生掌握利用计算机进行科学计算的基本理论和基本方法，并且学会将基本理论和基本方法应用于软件开发以及软件研制。2、本课程基本要求（1）掌握方法的基本原理和思想。（2）掌握方法处理的技巧及与计算机的结合。（3）掌握误差分析，收敛性及稳定性的基本理论。（4）学会进行可靠的理论分析，对近似计算要确保精度要求，要进行误差分析。（5）通过例子，学习使用各种计算方法解决实际计算问题。（6）通过上机实践，能编写算法和实现算法。（7）掌握数值计算中一些最基本、最常用的计算方法和算法。3、本课程与各课程的关系：由于本课内容包括了微积分、代数、常微分方程的数值方法，学生必须掌握这几门课的基本内容才能学好这一课程，同时，学习此课程还必须具备计算机系统的初步知识，掌握一门常用的高级语言，如：BASIC、PASCAL、C语言等，并须具备一定的编程能力。4、本课程的特点：（1）面向计算机，要根据计算机特点提供实际可行的有效算法。即算法只能包括加、减、乘、除运算和逻辑运算，是计算机能直接处理的。（2）有可靠的理论分析，能任意逼近并达到精度要求，对近似算法要保证收敛性和数值稳定性，还要对误差进行分析，而且都是建立在相应数学理论基础上的。（3）有好的计算复杂性。时间复杂性好是指节省时间；空间复杂性好是指节省存储量。这也是建立算法时要研究的问题，因为它关系到算法能否在计算机上完成。（4）要有数值实验。即任何一种算法除了从理论上要满足上述三点外，还要通过数值实验证明是行之有效的。计算方法最基本的立足点是容许误差，在误差容许的范围内对某一数学问题进行近似计算，得到能满足要求的近似结果。现实世界中误差是普遍存在的，由于世界上没有绝对精确的量具（绝对精确的量具是没有刻度的），因此人类通过量具采集的数据都是近似值，另一方面，我们的生产、实验工具都不是绝对精确的，这就使得人类在生产和科学实验中必需容许误差。计算机的应用可以分为二个方面，即数值计算和非数值计算。利用计算机进行数值计算的过程如下图所示：SKIPIF1<0在上图中，计算方法的任务是：由建立的数学模型给出可编程并由计算机能完成的计算方法，然后编程和上机求解。由于计算方法是编程后可由计算机求解的近似计算方法，如何确保近似解的精度显得尤为重要，必须深入讨论有关误差的基本概念和基本理论，为近似计算的精度分析打下基础。1、误差的来源（种类）误差的来源主要有以下四种（1）模型误差：建立数学模型时的误差。例如：在求重量的数学模型G=m*g中，重量G不是仅与质量和重力加速度有关，它还与温度、测量地点的海拔、地层结构等众多因素有关，为了使模型较为简单和实用，采用抓住主要矛盾的方法，去掉了大量对重量影响不大的次要因素，建立了上述重量的近似模型，由此产生了模型误差。（2）观测误差：采集数据时的误差。采集数据时，通常是依靠仪器和量具，由于没有绝对精确的仪器和量具，因此采集的数据有误差，此误差称为观测误差。（3）舍入误差：由于计算机字长有限而产生的误差。硬件再发展，计算机的字长总是有限的，在计算过程中，当数据的长度超过了计算机的字长时，计算机就会进行四舍五入，由此产生的误差称为舍入误差。（4）截断误差：无限形式的有限化而产生的误差。在计算中有时会运用无限形式的计算公式，例如台劳公式：SKIPIF1<0显然此公式无法进行计算，因此必需根据实际需要，从某一项起将后面的各项截断，即SKIPIF1<0由此产生的误差称为截断误差。§1.2绝对误差与相对误差、有效数字为描述方便，首先约定x*是精确值x的近似值。引入误差的概念，其目的是为了衡量近似值x*的好坏。（1）绝对误差：x*x由于精确值x通常无法确定，因此绝对误差无法计算，由此引入绝对误差限的概念。绝对误差限：绝对误差的一个上界。即：若|x*x|e，则称e为x*的绝对误差限。绝对误差限的性质是：A.不唯一这是因为|x*x|的上界是不唯一的。B.可确定只要我们对x*的实际背景有一定的了解，就不难确定|x*x|的上界。例如，x*表示身高，则|x*x|的上界可为3米。当x*是你求出的，那么为了说明你的工作认真，你一定会将|x*x|的上界估计得尽量小，因此在这种意义上绝对误差限可用来衡量x*的好坏。由于绝对误差限没有考虑问题的规模，因此有时它也不能衡量x*的好坏。例如：x是地球与太阳的距离，y是分子中二个原子间的距离，若|x*x|1公里，|y*y|1厘米，则并不能说y*比（2）相对误差：(x*x)/x*相对误差限：相对误差绝对值的一个上界。3、有效数字这里我们必须搞清楚什么是有效数字以及如何确定x*有几位有效数字。（1）有效数字的定义若|x*-x|<x*的某一位的半个单位，则称x*精确到这一位，并从这一位开始，一直到前面第一个不为零的数都是x*的有效数字。此定义实际上定义了什么叫精确到某一位和什么叫有效数字。例如：若x*精确到小数点后第3位，即指|x*x|0.510-3。（2）有效数字的判定方法方法一：四舍五入此方法首先确定x*是由x的哪一位四舍五入产生的，然后从这一位的前一位开始一直到前面第一个不为零的数都是x*的有效数字。例1若x=0.872596,x*=0.87，求x*的有效位数。解：SKIPIF1<0x*是由x的小数点后第三位四舍五入产生的，所以x*有二位有效数字。注意，方法一判定有效数字很简单，但有时会失效。例如，若x=0.272987x*=0.273102，此时无法用方法一确定x*的有效位数，原因是x*不是由x四舍五入产生的，在这种情况下，必须用有效数字的定义来确定x*的有效位数。即方法二：用定义此方法首先计算|x*x|，再判断它小于等于x*的哪一位的半个单位，然后从近一位开始，一直到第一个不为零的数都是有效数字。例2若x=0.62073，x*=0.6207，确定x*的有效位数。解：因为|x*x|0.00030.5104，x*精确到小数点后第4位，所以x*有四位有效数字。例3若x=0.080199，x*=0.802，确定x*的有效位数。解：因为|x*x|=0.000010.5105，所以SKIPIF1<00.5103，推出x*有三位有效数字。例4若x=6.28936，x*=7.3132，确定x*的有效位数。解：|x*x|=0.023570.5101，所以x*有二位有效数字。§1.3近似数的简单算术运算§1.4数值计算中误差分析的一些原则为保证计算结果的高精度，在进行数值计算时应遵循下述几个原则。（1）在进行除法时，要避免除数的绝对值<<被除数的绝对值。①为什么要“避免”？若不“避免”，则除出的结果很大，由于计算机字长有限，它装不下，因此会进行四舍五入，一个很大的数进行四舍五入时舍去的部分也会很大，这会使舍入误差变大。②怎样“避免”？因为用户只关心最后的计算结果，当中间计算过程中出现了除数的绝对值<<被除数的绝对值时，就应该换一种计算方法，以避免这种情况的发生，以后我们将会针对具体的计算问题来讨论“避免”的方法。（2）在进行减法时，要避免二个相近的数相减。①为什么要“避免”？若不“避免”，就可能失去大量的有效数字，例如：若a=30001和b=30000都有五位有效数字，因为a-b=1，所以结果至多有1位有效数字。②怎么“避免”？“避免”的思路与第1个原则中“避免”的思路相同，须针对具体计算问题来讨论。（3）要防止“大数吃小数”①什么是“大数吃小数”？我们用一个例子为说明。计算8756294874SKIPIF1<0，其中n=1020，0<ai<106。此题是一个很大的数与很多很小的数相加，若采用将大数依次与a1,a2,,an相加，由于计算机字长有限，因此在与ai相加时会进行四舍五入将ai舍去，这样，最后的结果仍是大数，这就是大数将a1,a2,,an吃掉了。②为什么要“避免”？尽管每个小数都很小，但它们很多，可能它们的和比大数还大，而最后计算工结果为大数，显然误差可能很大。③怎样“避免”？有的同学提出先将小数相加，然后再与大数相加，这个思路是对的，但有一个漏洞，因为小数相加到一定程度也会变成大数，它也开始吃小数了。可以采取分部相加的方法解决。第2章非线性方程（组）的近似解法§2.1引言方程f(x)=0的解称为方程的根。也叫做函数f(x)的零点。方程求根大致包括三个问题（1）方程有没有根？如果有根，有几个根？（2）哪里有根？求有根的区间，区间内的任意一点作为根的近似值。（3）根的精确化，已知一个根的近似值后设法逐步把根精确化，直到足够精确为止。本课程主要研究问题（2）和（3）。§2.2根的隔离求方程f(x)=0的解的近似值时，首先要确定若干个区间，使每个区间内只有的一个根，这个步骤称为根的隔离。对一般的方程，根的隔离有两种方法（1）试值法。求出f(x)在若干点上的函数值，观察函数值符号变化的情况，从而确定隔根区间。（2）作图法。画出y=f(x)的草图，观察曲线y=f(x)与x轴交点的大致位置，从而确定隔根区间。例1.2.1讨论方程f(x)=2x3-4x2+4x+2=0-1012-10-5-1012-10-5051011.52-1-0.500.5§2.3对分法11.52-1-0.500.5设有方程f(x)=0在（ab）内有且仅有一个根x*，这时有f(a)f(b)<0可用对分法求x*的近似值，方法如下（1）准备：计算区间（ab）两个端点的函数值f(a)，f(b)（2）对分：取c=(a+b)/2为（ab）的中点，计算f(c)

（3）判断：如果f(c)=0，则c为f(x)=0的根，否则检验： 若f(c)f(a)<0，则方程的根位于[ac]内，用c代替b，若f(c)f(b)<0，则方程的根位于[cb]内，用c代替a。

（4）检验：若|b-a|<e(e为精度要求)此时计算结束x*=c，否则转（2）。例1.3.1用对分法求方程f(x)=x3+2x-5=0在[12]内的根，[e=10-5有根区间f=inline('x^3+2*x-5')f(1),f(2)f=inline('x^3+2*x-5')f(1),f(2)fplot(f,[12]),gridon1.00001.50001.25001.50001.25001.37501.31251.37501.31251.34381.32811.34381.32811.33591.32811.33201.32811.33011.32811.3291方程的解x=1.3286§2.4迭代法设有方程f(x)=0在[ab]上有且仅有一个根x*，可用迭代法求x*的近似值，方法如下（1）将方程f(x)=0写成迭代形式x=(x)

（2）在[ab]上任取一个初始值x0。

（3）计算x1=(x0)

（4）若|x1x0|<e(e为精度要求)，此时计算结束x*=x1，否则令x0=x1转（3）。例1.4.1用迭代法解方程x=10x-2，x0=1分别采用迭代格式x=10x-2和x=log(x+2)，观察两个计算过程的区别。e=1e-3迭代过程:1.00000.47710.39390.37910.37640.3759迭代6次x=0.3759f=inline('f=inline('log10(x+2)')x=1x=f(x)例1.4.2用迭代法求方程f(x)=x3+2x-5=0的根，x0=1[SKIPIF1<0]。迭代过程:1.00001.44221.28371.34491.32201.33061.32741.32861.3281f=inline('(5-2*x)^(1/3)')x=1x=f(x)迭代9次x=1.3281§2.5牛顿迭代法f=inline('(5-2*x)^(1/3)')x=1x=f(x)牛顿法是解方程f(x)=0的重要方法，它也是一种迭代法。设有方程f(x)=0在[ab]上有且仅有一个根x*，可用牛顿法求x*的近似值，方法如下（1）求函数f(x)的导函数f(x)，牛顿法迭代公式为x=xf(x)/f(x)

（2）在[ab]上任取一个初始值x0。

（3）计算x1=x0f(x0)/f(x0)

（4）若|x1x0|<e(e为精度要求)，此时计算结束x*=x1，否则令x0=x例1.5.1用牛顿法解方程f(x)=x3-2x2-4x-7=0在[34]内的根[x0=4迭代过程:4.00003.67863.63293.6320迭代4次x=3.6320f=inline('x^3-2*x^2-4*x-7')fd=inline('3*x^2-4*x-4')x=4x=x-f(x)/fd(x)§2.6弦截法f=inline('x^3-2*x^2-4*x-7')fd=inline('3*x^2-4*x-4')x=4x=x-f(x)/fd(x)弦截法也是一种是解方程f(x)=0的迭代法，它的特点是不需要计算f(x)的函数f(x)，且收敛速度也相当快，是工程计算中常用的算法之一。设有方程f(x)=0在[ab]上有且仅有一个根x*，可用弦截法求x*的近似值，方法如下（1）求函数f(x)在区间[ab]的两个端点的函数值f(x0),f(x1)，其中a=x0,b=x1（2）计算x2=x1f(x1)[x1x0]/[f(x1)f(x0)]

（3）若|x2x1|<e(e为精度要求)，此时计算结束x*=x2，否则令x0=x1x1=x2转（§1.7用牛顿法解方程组设有非线性方程组u(x,y)=0，v(x,y)=0，在(xn,yn)按台劳级数展开，取展开式的第1，2项得到SKIPIF1<0SKIPIF1<0其中(xn,yn)是根的第n次近似值，如果Jn0方程组的第n+1次近似值(xn+1,yn+1)可用以下公式计算SKIPIF1<0例1.7.1用牛顿法解方程组u(x,y)=x3+y34=0，v(x,y)=x4+y23=0迭代初值x0=1,y0=1.4。（u/x=3x2u/y=3y2v/x=4x3v/y=2y例1.7.2用牛顿法解方程组u(x,y)=2x3y21=0，v(x,y)=xy3y4=0迭代初值x0=1.2,y0=1.7。（u/x=6x2u/y=2x v/x=6y3v/y=2xy21）例1.7.3用牛顿法解方程组u(x,y)=xcos(y)=0，v(x,y)=ysin(x)=0迭代初值x0=0,y0=0。（u/x=1 u/y=sin(y)v/x=cos(x) v/y=1）本章小结为了比较各种迭代方法的收敛速度，我们引入收敛阶的概念。设迭代过程xn+1=(xn)收敛于方程x=(x)的根x*，令en=xnx*，en称为迭代误差，如果存在实数P1和非零常数K，使得SKIPIF1<0，则称该迭代过程为P阶收敛的。P=1称为线性收敛，P>1称为超线性收敛，P=2称为平方收敛，显然P越大，迭代过程收敛的越快。可以证明当x*是方程f(x)=0的单根时，牛顿法是平方收敛的。当x*是方程f(x)=0的重根时，牛顿法仅为线性收敛。弦截法的收敛阶P=1.618。对分法的收敛速度与公比为1/2的等比级数相同。牛顿法：收敛速度最快，但要计算f(x)的导函数，计算量大，有发散问题。弦截法：收敛速度次之，不需要计算f(x)的导函数计算量比牛顿法小，有发散问题。对分法：收敛速度最慢，但简单有效，不存在发散问题。它一定收敛到有根区间[ab]内的某个根。第3章线性方程组的解法§3.1引言在科学实验和工程设计中，经常用到解线性方程组的问题。本章讨论用计算机求解线性方程组的两类主要方法：直接法和迭代法。解线性方程组的一般表达式SKIPIF1<0根据矩阵的性质可以写成SKIPIF1<0

简记为Ax=b其中SKIPIF1<0方程组Ax=b有唯一解的充分必要条件是|A|0。我们只讨论这种情况下的解法。解线性方程组的方法可以分为两类：一类是直接法，它只包含有限次的四则运算，在每次运算都无舍入误差的情况下，所得到的是方程组的准确解。由于实际计算中总是有舍入误差，所以实际得到的也是近似解。令一类是迭代法，它首先选择一组初始值，再运用同样的计算步骤，重复计算，得到近似解。由于这类方法中出现了极限过程，必须研究迭代过程的收敛性。本章主要介绍：直接法中的高斯消去法和主元高斯消去法。迭代法中的简单迭代法和塞德尔单迭代法。§3.2高斯消去法以n=4为例说明高斯消去法的计算过程，设有线性方程组SKIPIF1<0SKIPIF1<0SKIPIF1<0SKIPIF1<0经过3次消元步骤，得到以上形式。从最后一个方程中解出x4,依此回代得到方程组的全部解。6x1+3x2+2x3=6

例2.2.3用高斯消去法解方程组10x1+5x2+6x3=0

8x1+5x2+3x3=0方程组的增广矩阵[A|b]6326105608530消元6.00003.00002.00006.0000002.6667-10.000001.00000.3333-8.0000方程组系数矩阵主对角线元素为零，消元过程无法进行!例2.2.4用列主元高斯消去法解例2.2.3中的方程组。方程组的增广矩阵[A|b]6326105608530选主元1056063268530消元10.00005.00006.0000000-1.60006.000001.0000-1.80000选主元10.00005.00006.0000001.0000-1.8000000-1.60006.0000消元10.00005.00006.0000001.0000-1.8000000-1.60006.0000回代得到方程组的解5.6250-6.7500-3.7500§3.3矩阵的LU分解§3.4对称矩阵的LDLT分解§3.5线性方程组解的可靠性§3.6简单迭代法设有方程组Ax=b，变为迭代形式，x=Mx+f，或x(k+1)=Mx(k)+f，任取初始值x(0)程迭代得到x(0)，x(1)，x(2)，，x(k)，若极限SKIPIF1<0存在，则x*就是原方程组的解。以n=4为例SKIPIF1<0x(k+1) M x(k)f写成分量形式SKIPIF1<0定理1若SKIPIF1<0，则简单迭代法对任意初始值x(0)和f都收敛。定理2若SKIPIF1<0，则简单迭代法对任意初始值x(0)和f都收敛。定理3迭代公式x(k+1)=Mx(k)+f，对任意初始值x(0)和f都收敛的充分必要条件是矩阵M的各个特征值的模都小于1。§2.6雅可比迭代法与高斯-塞德尔迭代法在简单迭代法的基础上作改进x(k+1)=M1x(k+1)+M2x(k)+f，以n=4为例SKIPIF1<0x(k+1)M1x(k+1)M2x(k)f写成分量形式SKIPIF1<0定理1若SKIPIF1<0，则塞德尔迭代法对任意初始值x(0)和f都收敛。定理2若SKIPIF1<0，则塞德尔迭代法对任意初始值x(0)和f都收敛。定理3迭代公式x(k+1)=M1x(k+1)+M2x(k)+f，对任意初始值x(0)和f都收敛的充分必要条件是矩阵(I-M1)-1M松弛法（SuccessiveOverRelaxationMethod）x(k+1)=x(k)+(b-Ax(k))称为松弛因子，>1超松弛法，>1超松弛法，>1低松弛法。定理3松弛法对任意初始值x(0)和f都收敛的必要条件是0<<2。例2.6.1分别用雅可比迭代法和塞德尔迭代法解方程组SKIPIF1<0误差e<10-3雅可比迭代法迭代公式x(k+1)=Mx(k)+f,写成分量形式SKIPIF1<0初始值k=0(000),迭代过程x1(k) x2(k) x3(k) x1(k) x2(k) x3(k)-2.40005.00000.3000-4.00023.00311.9999-4.46124.24952.2802-4.00123.00002.0010-4.55582.74462.4671-4.00022.99922.0002-3.99132.62752.0345-3.99972.99981.9998-3.85792.98491.8865M=[0-0.4-0.2;0.250-0.5;-0.20.30]f=M=[0-0.4-0.2;0.250-0.5;-0.20.30]f=[-2.4;5;0.3]x=[0;0;0]x=M*x+f-4.03033.02372.0219-4.01382.98152.0132-3.99522.99001.9972-3.99543.00261.9960塞德尔迭代法迭代公式x(k+1)=M1x(k+1)+M2x(k)+f写成分量形式SKIPIF1<0clearM1=[0-0.4-0.2;00-0.5;000]M2=[000;0.2500;-0.20.30]f=[-2.4;5;0.3]x=[0clearM1=[0-0.4-0.2;00-0.5;000]M2=[000;0.2500;-0.20.30]f=[-2.4;5;0.3]x=[0;0;0]B=inv((eye(3)-M1))x=B*M2*x+B*fx1(k) x2(k) x3(k)-2.40005.00000.3000-4.40004.85000.3000-4.00403.07661.8667-3.99962.98712.0238-4.00003.00211.9961-4.00002.99972.0006-4.00003.00002.0000例2.6.2分别用雅可比迭代法和塞德尔迭代法解方程组SKIPIF1<0误差e<10-3（1）雅可比迭代法迭代公式x(k+1)=Mx(k)+fSKIPIF1<0取初始值k=0(000),迭代过程x1(k) x2(k) x3(k)1.00000-2.00005.0000-0.99600.0080-1.00805.00806.0080-1.00005.00006.0000-1.00005.00006.0000 分析，M的特征方程SKIPIF1<0（2）塞德尔迭代法迭代公式x(k+1)=M1x(k+1)+M2x(k)+fSKIPIF1<0，不收敛分析，B=(I-M1)-1M2的特征方程SKIPIF1<0例2.6.3分别用雅可比代法和塞德尔迭代法解方程组SKIPIF1<0误差e<10-3（1）雅可比迭代法迭代公式x(k+1)=Mx(k)+fSKIPIF1<0分析，M的特征方程SKIPIF1<0，所以简单迭代法不收敛。塞德尔迭代法迭代公式x(k+1)=M1x(k+1)+M2x(k)+fSKIPIF1<0分析，B=(I-M1)-1M2的特征方程SKIPIF1<0，收敛取初始值k=0(000),迭代过程x1(k) x2(k) x3(k)2.00001.000001.49800.2500-0.74802.24900.2495-1.43702.59380.4216-1.59392.58610.5039-1.54312.51960.5117-1.49902.49370.5027-1.49172.49450.4986-1.49682.49910.4988-1.50012.50060.4997-1.50062.50040.5001-1.5002本章小结本章讨论了解线性方程组的直接解法和迭代解法。直接解法比较适用与系数矩阵稠密（既零元素较少）的中、小型线性方程组，但对系数矩阵是带状或近似带状的大型线性方程组也适用。直接解法中的列主元高斯消去法具有精度较高和省时的优点，是计算机中常用的算法。迭代解法中主要介绍了雅可比迭代法、高斯-塞德尔迭代法和松弛法。迭代法具有计算公式简单、程序设计容易、占用计算机内存较少的优点。适用于解大型稀疏矩阵（既零元素较多）线性方程组。高斯-塞德尔迭代法是在雅可比迭代法的基础上改进得到，在很多情况下可以加快收敛速度，但它的收敛域与雅可比迭代法不同，因此不能互相取代。松弛法可以加速迭代过程的收敛速度，但要适当选择松弛因子（0<<2）。在选择迭代法时，要特别注意检验方法的收敛性问题。第4章矩阵特征值与特征向量的计算§4.1引言求矩阵的特征值和特征向量，是代数计算中的重要问题。在自然科学和工程中的许多问题，例如电磁振荡、桥梁的振动，机械振动等都可以归结为求矩阵的特征值和特征向量问题。矩阵A的特征值和特征向量是指，如果数和非零列向量x满足关系式Ax=x，则数称为A的特征值非零列向量x称为A的与特征值对应的特征向量。计算n阶矩阵A的特征值，就是求特征方程|AI|=0的根i(i=1,2,,n)。齐次线性方程组(AiI)x=0的非零解xi，是i对应的特征向量。本章讨论一些在计算机上计算矩阵的特征值和特征向量的较为稳定的数值算法。§4.2幂法和反幂法1、幂法：计算n阶矩阵A的模最大的特征值（主特征值）及对应的特征向量。任取n维列向量x(0)，用迭代公式x(k+1)=Ax(k)计算得到x(0)，x(1)，x(2)，设x(0)=a1v1+a2v2++anvn，因为Avj=ivj所以x(1)=Ax(0)=a11v1+a22v2++annvnx(2)=Ax(1)=a112v1+a222v2++ann2vn一般地有x(k+1)=Ax(k)=a11kv1+a22kv2++annkvn=1k[a1v1+a2(2/1)kv2++an(n/1)kvn]当k充分大时x(k+1)a11k+1v11x(k)向量x(k+1)与x(k)向近似地只差一个倍数，这个倍数就是模最大的特征值1。1、反幂法：Ax=x，A-1Ax=A-1x，A-1x=-1x，即A的特征值的倒数-1是A的逆矩阵A-1的特征值。用幂法求A-1的模最大的特征值，它的倒数就是A的模最小的特征值。§4.3雅可比方法对于2阶方阵SKIPIF1<0SKIPIF1<0令SKIPIF1<0对于n阶方阵（以n=3为例）SKIPIF1<0令SKIPIF1<0作变换矩阵SKIPIF1<0则有SKIPIF1<0令SKIPIF1<0作变换矩阵SKIPIF1<0则有SKIPIF1<0令SKIPIF1<0作变换矩阵SKIPIF1<0则有SKIPIF1<0一般地说，令SKIPIF1<0，可以将A中的元素的aij和aji变为0。在实际计算中采用以下公式SKIPIF1<0例4.3.1 用雅可比求对称矩阵SKIPIF1<0的特征值和特征向量。消去第i行第j列的元素[ij]=[12]A->1.00000.0000-0.70710.00003.0000-0.7071-0.7071-0.70712.0000消去第i行第j列的元素[ij]=[13]A->0.6340-0.32510.0000-0.32513.0000-0.62800.0000-0.62802.3660消去第i行第j列的元素[ij]=[12]A->0.59010.0000-0.08390.00003.0438-0.6223-0.0839-0.62232.3660消去第i行第j列的元素[ij]=[13]A->0.5862-0.02930.0000-0.02933.0438-0.62160.0000-0.62162.3700消去第i行第j列的元素[ij]=[23]A->0.5862-0.0252-0.0150-0.02523.41400.0000-0.01500.00001.9998消去第i行第j列的元素[ij]=[12]A->0.58590.0000-0.01500.00003.41420.0001-0.01500.00011.9998§4.4QR方法*QR方法是求一般矩阵A的全部特征值和特征向量的一种迭代方法。其基本思路是利用矩阵A的QR分解，通过迭代格式SKIPIF1<0将A化为相似的上三角矩阵（或分块上三角矩阵），从而求出A的全部特征值和特征向量。例3.4.1 用QR方法求矩阵SKIPIF1<0的特征值和特征向量。特征向量-0.8165-0.6215-0.6760-0.4082-0.6215-0.67600.40820.4770-0.2935特征值 1.00000.69724.3028本章小结本章介绍了求矩阵的特征值和对应的特征向量的几种方法。幂法可以求出矩阵的主特征值和对应的特征向量，优点是算法简单，但当|1/2|1时，收敛速度很慢。反幂法可以求出矩阵的模最小的特征值和对应的特征向量。雅可比方法是利用一系列正交相似变换（即平面旋转变换）把实对称矩阵A化为对角阵（近似），从而求出实对称矩阵全部特征值。QR方法是用镜向反射阵将矩阵A作QR分解，是一种求矩阵的全部特征值的有效方法。第5章插值与拟合§5.1引言已知表格函数y=f(x)xix0x1x2xn-1xnf(xi)y0y1y2yn-1yn构造一个公式p(x)近似地表示f(x)，解决这个问题的方法有两类：一类是插值法，另一类是拟合法，又称为逼近法。已知函数y=f(x)在互异点x0 ,x1,x2,,xn-1,xn上的函数值y0,y1,y2,,yn-1,yn，构造一个函数p(x)使得p(xi)=yi这样的问题称为插值问题。y=f(x)称为被插值函数，[x0xn]称为插值区间，p(x)称为插值函数，x0 ,x1,x2,,xn-1,xn称为插值点，在插值区间内部用p(x)代替f(x)称为内插，在插值区间外部用p(x)代替f(x)称为外推，R(x)=f(x)-p(x)称为插值函数p(x)的误差。§5.2插值多项式的存在和唯一性定理：给出n+1个插值点及函数值xix0x1x2xn-1xnf(xi)y0y1y2yn-1yn求一个n次多项式pn(x)=a0+a1x+a2x2++anxn(x0,y0)(x1,y1)(x2,y2)(x3,y(x0,y0)(x1,y1)(x2,y2)(x3,y3)(x4,y4)§5.3拉格朗日插值多项式1、给出2个插值点(x0,y0)，(x1,y1)可以得到一次多项式SKIPIF1<02、给出3个插值点(x0,y0)，(x1,y1)，(x2,y2)可以得到二次多项式SKIPIF1<0不难验证p2(x)满足插值条件p2(x0)=y0p2(x1)=y1p2(x2)=y23、给出n+1个插值点(x0,y0)，(x1,y1),,(xn,yn)可以得到一个n次多项式SKIPIF1<0其中SKIPIF1<0例5.3.1按下列表格求y(-0.5)和y(0.5)的值。x|123 y|7解：插值多项式SKIPIF1<0L2(-0.5)=0.2500L2(0.5)=4例4.3.2按下列表格求y(2.5)、y(4.5)、y(5.5)的值。x|解：插值多项式L4(x)=-0.7917x4+9.25x3-37.21x2+60.75x-32L4(2.5)=0.9180,L4(4.5)=6.1323,L4(5.5)=-8.9637§5.4均差插值公式已知函数f(x)在互异点x0,x1,,xn上的值为f(x0),f(x1),,f(xn)称SKIPIF1<0为函数f(x)在点xi,xj处的一阶均差。称SKIPIF1<0为函数f(x)在点xi,xj,xk处的二阶均差。一般地称SKIPIF1<0为函数f(x)在点x0,x1,,xn上的n阶均差。牛顿插值公式Nn(x)=f(x0)+(x-x0)f(x0,x1) +(x-x0)(x-x1)f(x0,x1,x2)+ +(x-x0)(x-x1)(x-xn-1)f(x0,x1,,xn)牛顿插值公式具有递推关系Nk+1(x)=Nk(x)+(x-x0)(x-x1)(x-xk)f(x0,x1,,xk+1)新增加一个节点，只需要增加计算一项(x-x0)(x-x1)(x-xk)f(x0,x1,,xk+1)§5.5差分等距结点插值公式1、差分已知函数f(x)在等距节点x0,x1,,xn上的值xx0x+hx+2hx+nhf(x)y0y1y2yn称△yk=yk+1yk为函数f(x)在点xk处的一阶差分。称△2yk=△yk+1△yk为函数f(x)在点xk处的二阶差分。一般地称△myk=△m-1yk+1△m-1yk为函数f(x)在点xk处的m阶差分。2、等距结点插值公式牛顿向前插值公式SKIPIF1<0牛顿向后插值公式SKIPIF1<0§5.6爱尔米特插值公式定理：给出n+1个插值点上的函数值和导数值xix0x1x2xn-1xnf(xi)y0y1y2yn-1ynf(xi)y0y1y2yn-1yn求一个2n+1次多项式H(x)满足Hn(xi)=yi，Hn(xi)=yi(i=0,1,2,,n)。SKIPIF1<0其中SKIPIF1<0SKIPIF1<0§5.7分段低次插值§5.8三次样条函数在xoy平面上给定n+1个点(x0,y0),(x1,y1),,(xn,yn)，构造一个函数S(x)满足以下条件（1）S(xi)=yi(i=0,1,2,,n)

（2）在区间(x0,xn)内S(x)具有连续的二阶导数（3）在每个子区间[xi1,xi]上S(x)（表达式是Si(x)）是一个三次多项式。满足以上条件S(x)称为三次样条插值多项式。三次样条插值多项式的推导设在子区间[xi1,xi]上S(x)=Si(x)(i=1,2,,n)由条件（1）得Si(xi-1)=yi1，Si(xi)=yi设S(x)在节点xi1处的二阶导数为Mi1，在节点xi处的二阶导数为Mi，即Si(xi-1)=Mi1，Si(xi)=Mi，显然Si(x)是x的线性函数，根据拉格朗日插值公式有SKIPIF1<0，记xi-xi1=hi有SKIPIF1<0将上式积分两次得到SKIPIF1<0利用Si(xi-1)=yi1，Si(xi)=yi定出积分常数C1和C2SKIPIF1<0解此方程组得到SKIPIF1<0代入上式整理后得到SKIPIF1<0SKIPIF1<0SKIPIF1<0类似地有SKIPIF1<0SKIPIF1<0因为Si+1(xi)=Si(xi)所以有SKIPIF1<0整理后得到关于位知数M0,M1,,Mn的线性方程组aiMi-1+2Mi+biMi+1=di(i=1,2,其中SKIPIF1<0边界条件：n-1个方程组，n+1个位知数，所有需要补充两个条件，常见的是以下两种（1）给定端点的二阶导数M0=a，Mn=b，特别地，当a=b=0时称为自然样条。（2）给定端点的一阶导数S(x0)=a，S(xn)=b，这时有SKIPIF1<0例5.7.1给定插值条件，和两种边界条件(1)m0=1,m3=0(2)M0=1,M3=0

x|0123

y|0000

m0=1,m3=0s1(x)=0.73333x3-1.7333x2+x0<x<1s2(x)=-0.2x3+1.0667x2-1.8x+0.933331<x<2s3(x)=0.066667x3-0.53333x(2)M0=1,M3=0s1(x)=-0.21111x3+0.5x2-0.28889x0<x<1s2(x)=0.055556x3-0.3x2+0.51111x-0.26667§5.9曲线拟合的最小二乘法给出表格函数xix0x1x2xn-1xnf(xi)y0y1y2yn-1yn求一个m(m<n)次多项式为pm(x)=a0+a1x+a2x2++amxm逼近这组数据。pm(x)的系数的确定：SKIPIF1<0即SKIPIF1<0(i=0,1,2,,n)SKIPIF1<0由微分学知，使(a0,a1,,am)达到极小值的a0,a1,,am满足必要条件SKIPIF1<0(K=0,1,2,,m)写成分量形式SKIPIF1<0正规方程组解此方程组得到a0,a1,,am即得到所要的多项式。例5.8.1有数据表如下，分别用二次和三次多项式逼近这组数据。x|-3-2-10123 y|p1(x)=0.1786x+0.5714p2(x)=0.17857x2+0.17857x-0.14286p3(x)=3.2895*10-17x3+0.17857x2+0.17857x-0.14286一次多项式逼近误差平方和2.8214，二次多项式逼近误差平方和0.1429，三次多项式逼近误差平方和0.1429。本章小结本章介绍了两种构造函数f(x)的逼近函数的方法—插值法和曲线拟合。关于插值法，讨论的是多项式插值，它要求所构造的多项式函数严格地通过给定的所有数据点。由Lagrange插值基函数的讨论，导出了Lagrange插值公式及其余项公式。作为对Lagrange插值的改进，建立了具有递推性的牛顿插值公式。Hermite插值是一种在插值节点上插值函数与被插值函数不仅有相同的函数值，而且还有相同的一阶导数值的多项式插值。由于高次插值会产生龙格现象，因此讨论了分段插值法，重点介绍了三次样条插值公式。关于曲线拟合，讨论的是曲线拟合的最小二乘法，它不要求所构造的逼近函数严格地通过给定的所有数据点，只是在多项式中，以使得残差的平方和最小为标准，选择多项式，以作为被逼近函数的近似替代。本章最后介绍了数值微分，其中5点公式是常用的。第6章数值积分和数值微分§6.1引言计算定积分的牛顿-莱布尼兹公式SKIPIF1<0，其中F(x)是函数f(x)的原函数。在实际应用中（1）常遇到某些函数f(x)的原函数不能用初等函数表示，例如sin(x2)，cos(x2),sin(x)/x，1/log(x)等。（2）有些函数f(x)是用表格形式表示，无法得到它们的原函数。（3）有些函数f(x)的原函数十分复杂，不利于工程上的使用。因此，要研究计算定积分的近似方法：数值积分。§6.2牛顿一科特斯型积分公式由上一章知，任意函数f(x)可以用一个拉格朗日插值多项式n(x)近似表示，因此f(x)的定积分也可以用n(x)的定积分近似表示SKIPIF1<0以此为基础得到的数值积分公式称为牛顿一科特斯型积分公式。SKIPIF1<0其中SKIPIF1<0取a=x0<x1<<xn=b为一组等距节点，令x=a+bt，h=(b-a)/n得到SKIPIF1<0记SKIPIF1<0，于是SKIPIF1<0所以得SKIPIF1<0对不同的n系数如下n牛顿一科特斯型积分公式系数余项（误差）11/2，1/2 (梯形公式)(1/12)h3f(2)()a<21/6，4/6，1/6 （抛物线公式或辛卜生公式）(1/90)h5f(4)(31/8，3/8，3/8，1/8(3/80)h5f(4)(47/90，32/90，12/90，32/90，7/90 （柯特斯公式）(8/945)h7f(6)(519/288，75/288，50/288，50/288，75/288，19/288(275/12096)h7f(6)(641/840，216/840，27/840，272/840，27/840，216/840，41/840(9/1400)h3f(2)(§6.3复合求积公式1、复合梯形公式取a=x0<x1<<xn=b为一组等距节点将积分区间[ab]分为n个子区间[xkxk-1]xk–xk-1=(b-a)/n=h，在每个子区间[xkxk-1]上应用梯形公式得：SKIPIF1<0于是得到复合梯形公式SKIPIF1<02、复合辛卜生公式在积分区间[ab]取点a=x0<x1<<x2n-1<x2n=b将[ab]分为2n个子区间，令x2k-2x2k=(b-a)/n=2h，在每个子区间[x2k-2x2k]上应用辛卜生公式得：SKIPIF1<0（k=1,2,,n）于是得到复合辛卜生公式SKIPIF1<03、复合柯特斯公式在积分区间[ab]取点a=x0<x1<<x4n-1<x4n=b将[ab]分为4n个子区间，令x4k-4x4k=(b-a)/n=4h，在每个子区间[x4k-4x4k]上应用柯特斯公式得：SKIPIF1<0（k=1,2,,n）得到复合柯特斯公式SKIPIF1<0（k=1,2,,n）4、步长h的自动选择复合梯形公式|T2nTn|<3e（e表示允许误差）其中Tn和T2n分别表示取n和2n时用复合梯形公式计算得到的积分近似值。（2）复合辛卜生公式|S2nSn|<15e（e表示允许误差）其中Sn和S2n分别表示取n和2n时用复合辛卜生公式计算得到的积分近似值。（3）复合柯特斯公式|C2nCn|<63e（e表示允许误差）其中Cn和C2n分别表示取n和2n时用复合柯特斯公式计算得到的积分近似值。§6.4龙贝格求积公式梯形公式与辛卜生公式之间的关系SKIPIF1<02、辛卜生公式与柯特斯公式之间的关系SKIPIF1<03、龙贝格积分公式SKIPIF1<0§6.5高斯求积公式一个求积公式，若对于任何次数不超过m的多项式都准确成立，则称这个求积公式的代数精确度为m。定义：使插值求积公式SKIPIF1<0的代数精确度为2n-1的节点x1,x2,,xn称为高斯点，对应的插值求积公式称为高斯求积公式。定理1：若x1,x2,,xn是高斯点，则SKIPIF1<0定理2：高斯求积公式的系数恒为正，且有SKIPIF1<0n节点xk(n)系数Ak(n)余项（误差）积分区间[-11]102f(2)()/31<<120.5773503+0.577350311f(4)()/13530.5773503

0+0.57735035/98/9

5/9f(6)()/1575040.8611363

0.3399810

+0.3399810+0.86113630.3478548

0.6521452

0.6521452

0.3478548f(8)()/347287550.90617990.5384693

0

+0.5384693+0.90617990.2369269

0.4786287

0.5688889

0.4786287

0.2369269f(10)()/1237732650§6.6二重积分的数值积分法§6.7数值微分给出表格函数如下表所示，求f(x)在节点f(xi)处的导数的近似值。xix0x1x2xn-1xnf(xi)y0y1y2yn-1yn方法：用插值多项式的导数近似f(x)的导数。三点公式：SKIPIF1<0SKIPIF1<0三个相邻节点的选法，一般是在所考察的节点两侧各选取一个节点，如果一侧的无节点，则用另一侧的节点补足。一阶导数的