数值分析讲稿

1、数值方法（A）讲稿周 富 照数值方法（A）讲稿（周富照）数值方法研究的对象：各种数学问题求解的数值计算方法。就是研究用计算机解决数学问题的数值计算方法及其理论。包括方法的稳定性、收敛性及误差分析，还要根据计算机特点研究计算时间最省的计算方法。实际问题的解决过程：实际问题数学问题计算方法上机实现实际问题第一章 绪 论掌握误差的来源、衡量误差大小的一些标准、控制误差的一些原则1.1 误差的来源与分类 在工程和科学计算中，估计计算结果的精确度是很重要的。而影响精确度的是各种误差。误差按照它们的来源可分为以下四种：1. 模型误差数学模型一般只能近似的描述实际问题，如匀速运动认为速度(近似地)为常数，匀

2、加速运动认为速度为时间t的一次函数。在动力学中，若认为阻力R是速度 v = x (t) 的一次函数：R = v，将v的高次项的影响略去，可得线性微分方程模型(见同济高数)。2. 观测误差数学模型中包含的某些常数(时间，长度等)由于受仪器的限制观测结果不能绝对准确。如测量距离时为了减少测量误差，一般取两次往返测量的平均值。3. 截断误差有些函数需要用无穷级数计算，计算时只能取前几项。如 当x充分小时取 sin x x，其误差 。4. 舍入误差在计算机上要把一些数字四舍五入后再计算如：因为我们主要是在已给数学模型的基础上研究计算方法的，所以只考虑后两种误差。1.2 误差概念或误差的衡量标准绝对误差

3、（限）、相对误差（限）和有效数字 1. 绝对误差 设某量准确值为x，近似值为x*，称 e = x x*为近似值x*的绝对误差。简称为误差。 在同一量的不同近似值中，|e (x)|越小，x*的精确度越高。2. 绝对误差限 由于精确值x实际上是不知道的，故绝对误差e也不能求出。在实际计算中，可根据情况事先估计出它的大小范围。即预先指定适当小的正数，使得 | e | = | x x* | 称为近似值x*的绝对误差限。有时也用 x = x* 表示近似值的精确度或准确值的范围。 例如，取的近似值为 a，则 e = a， | e | 0.0003。 3. 相对误差绝对误差有时不足以表示近似值的好坏。例如，

4、若有两个近似数x1 =100 1（m）， x2 =1000 1(m)其绝对误差限都是1(m), x1* = 100(m), x2* = 1000(m)，与近似值本身比较，x2*较精确。记 称er为近似值x*的相对误差。如以上问题中的x1相对误差为1，x2为。 相对误差越小者越精确。相对误差在误差分析中更能反映出误差的特征。它无量纲，与所取单位无关。4. 相对误差限 和绝对误差一样，er 也不能求出，可预先指定适当小的正数 r，使得 ，或 。 r称为近似值x*的相对误差限。5. 有效数字近似值的准确程度可用有效数字来描述。定义 设x*是的近似值，若x*的绝对误差限是它的某一位的半个单位，并且从该

5、位到它的第一位非零数字共有位，则用x*近似是具有位有效数字。相对误差与有效数字位数有关。有效数字位数越多，相对误差限越小。例1 设，0.005 是经四舍五入得到的近似数，分别求出它们的绝对误差限和相对误差限。由此可得什么结论? 解 绝对误差限：，。 相对误差限：0.02%，0.02%，10%。例2 下列各近似值的绝对误差度是，它们有几位有效数字。 （4位），（2位），（0位） 1.3 误差的传播设 是可微的,是的近似值, ，则问题是：当()很小时，是否很小？ 当()很小时，是否很小？ 下列公式很容易从对应公式求得：， ，故可取 ， 简单地说，两数之和(差)结果保留的小数位数与位数较少的加数相同

6、，两数之积(商)保留的数字位数与有效数字较少的因数相同。但其误差可能超过末位半个单位，这样保留的数除最后一位数字外都是有效数字。详细分析如下例. 例2 计算 ( 62.33)并分析误差，设因数都是四舍五入所得(计算结果只保留有效数字)。解 设x，y，由 d(xy) = ydx + xdy 得 (xy) = |y|(x) + |x|(y)。 (x，(y， (xy) + 0.0005 = 0.0016 + 所以要四舍五入保留1位小数(为什么?)： ( 62.33) = 20.1949 。1.4 数值运算中的一些原则 算法稳定性的若干原则一个算法若计算结果受原始数据误差影响很小，便称之有较好的数值稳

7、定性。因为在计算机中进行运算会产生舍入误差。为了减少舍入误差的影响，设计算法时应遵循以下一些原则：1. 要设法控制误差的传播2. 相近的数相减的改进措施因为 | x y|太小会使相对误差扩大。1)若是两个根式相减，可将分子分母同乘以它的有理化因式。如取四位数字计算： (只有一位有效数字) (有四位有效数字)2)若是两个三角函数相减，可和差化积后再计算。如：o o sin 30 oo o。3) 若是超越函数在两个相近点的函数值相减，可用Tailor公式或中值定理： f (x2) f (x1) = f ( ) (x2 x1)， (x1, x2)例 比较三种算法：ln 3.000 001 ln 3

8、= ln (3.000 001 / 3) 0.000 001 / 3.000 000 5。用计算器计算结果分别为: 3.333 3107，3.333 299 444107，3.333 332 778107 (把计算结果乘以10 7才能显示10位)。更精确的解为3.333 332 777 777 9107。3. 两数相加要防止较小的数加不到较大的数中连加时，要把绝对值较小的数放到前面先加。如：在4位有效数字的计算机上计算 0.3 + 0.3 + 0.4 + 1234，1234 + 0.4 + 0.3 + 0.3 结果不同。4. 减少运算次数这样可以减少误差积累。如计算多项式可用秦九韶算法： 例

9、2 x3 3 x2 + 4 x 5 = (2 x 3) x + 4) x 5。5. 避免除数绝对值太小因为 d (x/y) = (ydx xdy) / y2，| y | 太小会使绝对误差扩大。6. 编程时避免用等式条件 例题下列各题右端算法精度较高：1. ， |x| 1，和差化积。2.3.，N充分大。4.，|x| 1，用Tailor公式。记住 (1+x)a, ex, ln(1+x), sin x, cos x, tan x, cot x 的Tailor公式：，。 5计算，取，利用下列等式计算，哪一个得到的结果最好？6序列满足递推关系）若（三位有效数字），计算到误差有多大？这个计算过程稳定吗？

10、7 给定 ，试用四位数学用表求的近似值。 甲方法：， 乙方法：，而，问题那一个方案较正确？分析：记 ，其中 ，三角函数表给出了四位数字，它准确到小数后第三位，而第四位是经过“四舍五入”得到的，即有，乙比甲好，实际答案。第二章 插值与逼近2.1 引 言1、插值问题：在实践中经常遇到这样的函数，它太复杂或仅仅是一些测量数据，就需要用较简单函数近似表示: 。 已知y = C a，b 在n+1个互异节点a = x0 x1 xn = b处的函数值=，i = 0, 1, , n。为某一函数类，求满足: (xi) = yi （1）称这个问题为插值问题，为的插值函数，为被插值函数，x 0, x 1, , x

11、n 为插值结点，（1）式为插值条件。2、解决问题包含的内容（1）确定（有限维空间只要确定一组基）；（2）讨论满足（1）式的是否存在？是否唯一及如何求？（3）讨论截断误差 的表达式与估计。常用插值函数类有: 多项式，三角函数，有理函数等。对应插值法分别称为代数多项式插值，三角多项式插值，有理函数插值。下面主要讨论代数多项式插值。2.2 代数插值问题1、代数插值问题设表示所有次数不超过的多项式的集合。其它条件如上，求满足, i = 0, 1, , n （1）2、存在唯一性：定理 满足条件（1）式的代数插值多项式存在且唯一。，有 （2）求求一组数，将（2）式代入（1）式： （3）的存在唯一性方程组（

12、3）的解的，的存在唯一性（3）的系数矩阵的行列式为n+1阶Vandermonde行列式由于结点互异，所以 Vn 0，故方程组（3）存在唯一解(a 0，a 1，a n )。3、插值多项式的求法: Lagrange插值多项式下面我们用另一种方法确定多项式。令 ，k = 0, 1, , n （6）则lk (x)满足条件 k=0, 1, , n （7）记 （8）则有 ，i = 0, 1, , n.所以 称上式为Lagrange插值多项式, 称lk(x)为以、为结点的Lagrange插值基函数。 注：（1）n=1时，线性插值公式， ，或 （2）n=2时，二次插值（或抛物线插值）公式：4、插值余项设为 f

13、 (x) 的n次Lagrange插值多项式，其余项为R n (x) = f (x) 定理 设 f C n+1 a, b，则对任意 x (a, b) ，存在 (a, b) ，使得 （4）其中 。证 当x = xi 时定理成立。下设 x xi，作 （可根据不同层次进行取舍） （5）则 C n +1 a, b，且，即在a, b上有 n + 2个零点 t =。应用Rolle中值定理，在(a, b)内至少存在n + 1个点 0 , 1 , , n, 使得：F ( i) = 0，i = 0, 1, 2, , n。如此反复应用Rolle中值定理，可知在(a, b)内至少存在一点使得因为，所以有由此可得定理的

14、结论。 几点说明：(1) 插值多项式本身只与插值结点及f (x)在这些点上的函数值有关，而与函数f (x)并没有关系。但余项R n (x) 却与f (x)联系很紧。具体估计时f (n+1)()可取 |f (n+1)(x)| 在区间上的最大值：|Rn(x)| M |n+1 (x)| / (n+1)!，M = max x a, b |f (n+1) (x)|。有些函数的 |f (n+1)(x)| 随 n 增长很快，最好采用分段低次插值。 (2) 若f (x)为次数 deg (f ) n 的多项式，那么以n + 1个点为结点的插值多项式就一定是其本身。因此时 f (n +1) ( ) 0，故R n

15、(x) = 0。还进一步推出一些结果，见书后习题。(3) 当点x位于x 0 , x 1 , x n 的中部时， | n+1 (x) | 比较小，精度较高。位于两端时，精度较低。若点x位于 x 0 , x 1 , x n 的外部，一般称为外插，精度更低。见下图所示。 4 (x)= (x2 0.09) (x2 0.01) 6 (x)= (x2 0.25) (x2 0.09) (x2 0.01)图5-1图5-2例1 sin 15o，sin 16o，求 o解 故 o = 0.2622.例2 x和y的观察数据如下表，试求Lagrange插值多项式x012 y 123解 这说明的次数可以小于n 。例3 已

16、知函数的观测数据为x1234 y 0-5-63（1）试求拉格朗日插值多项式；（2）分别用线性插值和二次插值求解：（1）（a）求插值基函数 ， =（b）求(2) 线性插值：（如何取点？取） ，二次插值：（如何取点？取）以上结果说明什么？2.3 分段线性插值1是f(x)在a,b上的拉格朗日插值函数，当插值节点无限加密时，在a,b趋于f(x)吗？2分段线性插值 已知y = f (x) C a , b 在n+1个互异节点a = x0 x1 xn = b处的函数值yi = f (xi ), i = 0, 1, , n. 试求,，使 （1）， i = 0, 1, , n. （2）在，i = 0, 1, ,

17、 n-1上是一次多项式这样的称为f(x)在a,b上以a = x0 x1 xn = b为节点的分段线性插值多项式。注 （）这样的使存在唯一的：只将上的一次插值多项式拼接起来即可。，i = 0, 1, , n-1（）记，i = 0, 1, , n-1，定理 若f(x)在a,b上连续，则，。2.4 Newton插值公式 差商与Newton插值多项式Lagrange公式优点：对称，很容易写出。缺点：增加结点要重新计算。若写成 (1)则可以逐次计算出： a0 = f (x0), , 每增加一个结点只需增加一项，前面的项不变。下面利用差商概念导出Newton差商插值多项式，就能很快确定系数。 1. 差商设

18、函数y = f (x)在xi , xj 上有定义，称 为 f (x)在xi , xj 上的一阶差商。称 为 f (x)的二阶差商.已知 k 阶差商 f xi , xi+1, xi+k，f x i +1, x i+2 , x i+k+1 ，定义 k + 1 阶差商为 并规定零阶差商为函数本身：f xi = f (xi)2.5 埃米尔特（Hermite）插值1、问题：已知y = f (x) C a , b 在n+1个互异节点a = x0 x1 xn = b处的函数值yi = f (xi )和导数值, i = 0, 1, , n。求 ，使 ， i = 0, 1, , n （1）2、Hermite插值

19、多项式的求法 令 ， （2）其中是不超过次的多项式。显然 由 ，得，由多项式插值：代入（2）式得+。2.6 三次样条插值多项式插值问题：结点增加，次数升高，计算量增大，波动增大。图5-3为的图象(实线)及将1, 1分成10等分所作Lagrange插值多项式P10 (x)的图象(虚线)。P10 (x)只中部较好的靠近f (x)，越往两端震动越大(叫Ronge现象)。改进方法：分段低次插值。若还要求结点处要光滑，就是样条插值函数。1三次样条函数 设给定区间 a, b 上n + 1个点：a = x0 x1 xn = b。若 s (x) 满足：(1) 在每个子区间xk1, xk (k = 1, 2,

20、, n)上， s (x)为不超过三次的多项式；(2) 在a，b上有二阶连续可导则称s (x) 是a，b上以 x1, x2, , xn1 为内部结点的三次样条函数。a，b区间上以为节点的所有三次样条函数的全体记为。2三次样条插值问题 求 使，满足上式的s (x)称为的（在a，b上以为节点）三次样条插值函数。3三次样条插值函数的求法（之前简单分析三次样条插值函数的存在性和唯一性问题）（1）三弯矩方程组方法（1）令(= 0, 1, , n)，设x , ，则s (x)为过两点，的线性函数。设，则 （1）对上式积分两次，得 (2) (3)根据插值条件：，得 ，由此解得 代入(2)式得 (4)由此得x ,

21、 时有： (5)当x , 时， (6)根据三次样条插值函数的定义应有：，由此得： (7)令 , 代入(7)式得三弯矩方程，或 (8)此方程组有 n + 1个未知量，但只有 n 1个方程，还需根据实际问题补充两个条件(端点条件)。 已知两端点x0及xn处的导数和：，。由式(5)，(6) 也即 (9)其中 将(8)与(9)两个方程组联立，得 (10) 函数 y = f (x) 在端点曲率为已知常数：，=。代入方程组(8)即可。M0 = Mn = 0 时称所得的样条为自然样条，表示两端点为简支的情况。 函数 y = f (x) 以 b a = xn x 0为周期： y0 = yn , s(x 0)

22、= s(x n)， s(x 0) = s(x n)，由于 hn+1 = h1，Mn = M 0，yn = y0，利用(5)式和(6)式可得 。于是 由此得 (11)其中 由(8)式和(11)式得 (12)可用推广的追赶法求解。（2）三转角方程组方法设，由，利用Hermite插值公式有+， （1）由上式得类似得 由（）有 （） （2）考虑边界条件（1），即，已知，代入（2）式得 （3）（用什么方法求解？）解方程组（3）得 ，代入（1）式得。这样求的方法称为三转角方程组方法。若考虑边界条件（2），已知，得 (4)其中 ，结合（2）式和（4）式得 （5）解方程组（3）得 ，代入（1）式得。至于周期性

23、边界条件类似可得结论。例1 已知插值条件：，与边界条件：，求在1，3上的三次样条插值函数。解：（1）， 三转角方程组： 结合：，有（2）代入（1）式得 例2 给出四个样点(1, 1)，(2, 3)，(4, 4)，(5, 2)，求其各个子区间上的样条插值函数s(x)，并求f (3)。设已知 y0= 2，y3= 3。解 求Mk 的方程组为：列表计算,如下，其中Dk = (yk yk 1) / hk-1：khk= xk+1 xk x kykDk 0112 (=y0)012/31232 3 21/32441/2 5 3152 2 64 3 (=y3)代入方程组得 解得 代入(3)，化简得 而 f (3

24、) s (3) = (533 4 32 + 199 。例3 已知数据表x1245y1342及边界条件。求其各个子区间上的样条插值函数。解：（1）n=3，；，（2）三弯矩方程组为（3）代入(3)，化简得 。例4 已知插值条件：，与边界条件：。求在1,3上的三次样条插值函数。解 求的方程组为：列表计算，如下，其中Dk = (yk yk 1) / hk-1：khk= xk+1 xkx kykDk0111(=y0)611/21232-1821/213-1-430 31(=y2)代入方程组得 解得 ，代入(3)，化简得 2.7 最佳平方逼近与正交多项式一、几个概念1：上所有连续函数组成的集合2上的权函数

25、：（1），；（2）存在，（3）且非负，若，则=0，3上的内积：若，是上的权函数，则称为关于在上的内积。4若与的内积为零： 则称f与g关于在a,b上正交。5若满足则称为正交函数族。6，若当且仅当，时， 则称线性无关。7一般逼近问题 已知，求使其中表示在空间中与的距离。称为最佳逼近元素。若取 ，则上问题称为最佳一致逼近问题；若取，则上问题称为最佳平方逼近问题。二、最佳平方逼近问题1、最佳平方逼近问题解的存在唯一性设，且线性无关。， （1）求求一组数，记 （2）， （3）记 ，方程组（3）化为 （4）因线性无关，所以。定理 设如上定义，则对，存在唯一的使2、求最佳逼近元素的步骤：（1）令，求， 。（

26、2）求解方程组（4）得：（3）将代入（1）式即得。问题：（a）方程组（4）是否有解？如何求解？用什么方法求解？(b) 与有关，是否可选取使尽可能简单？2.8 正交多项式系1、定义 （1）Legendre多项式系 在-1,1上， ，（2）Chebyshev多项式系-1,1上，2性质（1）（正交性）Legendre多项式系是-1,1上以的正交多项式系Chebyshev多项式系是-1,1上以的正交多项式系证： （2）（奇偶性），（3）（递推公式），证： （4）（）在-1,1上有个相异的实零点。例1 设（2+x），x-1,1，求它在上的最佳平方逼近多项式。解：取，例2 设，求它在上的最佳平方逼近多项式

27、。 解： ， ， ， ，解方程组 ， 2.9 曲线拟合的最小二乘法1、问题 已知一组数据，求使 （1）求满足（1）式的的方法称为最小二乘法（曲线的最小二乘拟合），称为最小二乘拟合曲线。2、求解取 ，线性无关。， （2）记 求求一组数；求求的最小值点。由多元函数求最值的方法，令， （3）记 ，（2）式 （4）因为，线性无关，可以证明，从而（4）有唯一解，解记为 。可以证明总结求的步骤： （1）先确定，；（实际情况如何确定？） （2）计算 ， （3）解方程组（用什么方法求解？）得解 ； （4） 得 例1 已知 012312579421试求它的最小二乘拟合曲线（。解：（1）取直角坐标系，描点，由图可

28、知，这些点位于一条双曲线附近。取 ，即，； （2） ，=1.842857, =1.310408，=16， (3) 解方程组得解，问题：能否找到这样一组基使方程组（4）的系数矩阵为简单形式？可找到这样的使=（最小二乘的）多项式拟合若取，则总结求（最小二乘多项式拟合曲线）的步骤： （1）计算 ， （3）解方程组（用什么方法求解？）（称为正规方程）得解 ； （4） 得 若记 ，则 ，正规方程组 例1 用一次函数 y = a + b x 拟合右表所列数据123035解 这里 ，.(正规方程组：)解得 b = 5/2，a = 7/3。所求拟合曲线为 y = 2.5 x。例2 求一个经验函数 (x)=ae

29、bx 使它能与下表所列数据相拟合：xk12345678yk解 对经验公式两边取对数： ln (x) = ln a + bx。令A= ln a，B=b，u= ln (x)，则拟合曲线方程可化为：u=A+Bx。故只需对x和ln y的函数表用例1中的方法先求出A，B，再求a，b。其余请同学们自己完成。例3 已知三次样条函数在上的表达式为，且，求在及上的表达式。答案：例4 已知三次样条函数在上的表达式为，且，求在及上的表达式。答案：设在，及上的表达式分别记为，和；又设，由，得，。其余同理。例5 设是以互异的为节点的Lagrange插值基函数，证明 简答：（1）记 ，的Lagrange插值函数为，因，所

30、以由插值余项公式，即（2） 。例6 已知在，上，求。简答： ， ， ，。例7 已知插值条件：，与边界条件：。求在1,3上的三次样条插值函数。简答：（1）， ，； （2）三转角方程组 ，结合 ，得； （3）由公式+ +，算得。第四章 数值积分1 理论：若在上连续，且 。则2 实际问题：下列情况如何计算？ （1）的计算相当复杂时； （2）不是初等函数时； （3）不知的解析式，而只知在若干点上的值。3 机械求积法取：用形如 其中 只与有关，而与无关。这类方法称为机械求积法；：求积节点；：求积系数4.1 NewtonCotes求积法1NewtonCotes求积公式 设：，作在上拉格朗日插值多项式，记，

31、 （1）因为用近似，所以可用近似，于是得 ， （2）（2）式称为阶NewtonCotes共识，称C i(n) (i = 0, 1, , n)为Cotes求积系数。注： 满足 （可简单证明），； 当时，得梯形公式：当n=2时得抛物线(Sympson)公式：当n=3时， 例 用n=15阶 Newton-Cotes求积公式计算()计算结果n123452NewtonCotes求积公式的误差 截断误差， 代数精度定义 若某个求积公式对于任意次数不高于m的多项式准确成立，而存在一次数为 m + 1的多项式使之不准确成立，则称这一求积公式的代数精确度为m。定理 设f (x)在区间a, b上具有连续的二阶导数

32、，则梯形求积公式的误差为 设 f (x) C 4 a, b，则抛物线公式的误差为： n阶NewtonCotes求积公式至少有次代数精度。 当为偶数时，n阶NewtonCotes求积公式至少有次代数精度。证 这里n = 1，在区间a, b上2 (x)=(x a)(x b)不变号，f (x)连续，由积分第二中值定理，存在 (a, b)使得 因为NewtonCotes求积公式的余项为： 若deg (f ) n，则 f (n+1) ( ) = 0，故n阶NewtonCotes求积公式对于任何次数不高于n的多项式是精确的。 设为一次数为n+1的多项式，则（常数），所以 =0例 分别求梯形公式和抛物线公式

33、的代数精度（分别为1和3）。（注意求的步骤和方法）3NewtonCotes求积公式的数值稳定性与收敛性 当不大时，数值稳定； 当较大时，数值不稳定； 当（或）时，不收敛于。证： 设的实际计算值为，的近似值为，(b-a)， (b-a)， ，因，所以结论成立。4.2 复化求积公式为提高精度，使用高阶求积公式，这样存在的问题：（1）很大时，计算很复杂；（2）很大时，不一致收敛于，从而不收敛于；（3）很大时，数值不稳定。1复化梯形公式将a, b n等分， ，在每个子区间上应用梯形公式， （1） 得复化梯形公式：， （2） 复化梯形公式的截断误差： ， （3）（3）式的证明：由梯形公式的截断误差易知 ，

34、 （4）由连续函数界值定理知, (a, b) 使得： 代入（4）即得（3）式。实际估计取 2复化抛物线公式类似复化梯形公式，把区间a, b分成n个相等的子区间xi, xi+1(i = 0, 1, , n 1)，每个子区间上应用抛物线公式， （5） 复化抛物线公式：， （6） 复化抛物线公式的截断误差： ， （7）实际估计取 同理可推出复化Cotes公式。 例1 若用复化梯形公式和复化抛物线公式计算积分，问积分区间要等分多少才能保证有5位有效数字?解 由余项公式因为 f (x) = f (x) = = e x, b a = 1故 或 由 | R | 10 4/2，得 n (e10 4/6)1/2

35、 =68 (梯形公式)，或 n (e104/1440)1/4 =2.1(抛物线公式，因区间是2n，实际上还要乘以2，故6等分即可)。由此可知，使用抛物线公式比使用梯形公式可大大提高计算效率。例 分别用复化梯形公式和复化抛物线公式计算，使精确到解：（1）取，81689， （2）取，例2 将0, 1区间分成8等分，分别用复化梯形公式和复化抛物线公式计算：解 将分点及函数值列表如下i012345678xi011/(1+xi2)10.5 ，I ，I (1/24)1+0.5+2(0.9412+0.8+0.64)+4注：；复化梯形公式和复化Simpson公式都是数值稳定的。3 变步长公式 实际计算时因为难

36、以事先估计步长，往往借助于计算机自动选择步长。下面介绍变步长抛物线公式（抛物线法的递推化）。 第一步 将区间a, b分成2等分， 计算 第二步 将区间a, b分成4等分， 计算(xi = a + ih, h = (b a)/4，i = 0, 1, 2, 3, 4)： 若，则停止计算；否则继续将每个小区间2等分，计算，若， 则以作为所求积分近似值。否则继续将区间分半进行计算。注： （1）当在a,b上变化很小时，有 ； （8）（2）当在a,b上变化很小时，有 。 （9）（3）当和无法估计时，如何估计的误差？事后估计方法：用做为估计，即|（或）（或）。 由（8）和（9）式可以认为 比的效果好，比的效

37、果好。实际上可以推出， ，（Romberg求积法）（5），与的关系：，：， ， ， 。4.3 Romberg 求积法一、Richardson外推法设常量，用的函数来近似它，设已知+， （1）其中与无关，单调递增，取，+， （2） （3）有+， （4），（），与无关，一般地，可归纳做出， （5）以上这种方法称为Richardson外推法。二、Romberg求积法将Richardson外推法应用于复合梯形公式可得收敛性很好的数值求积法Romberg求积法。记 （6），当在a,b上满足一定条件时，有+， （7）其中与无关，（6）式中用代，+， （8） （9）有+， （10）类似从（9）式出发，+，

38、（11）一般地，可归纳做出， （12）由上式，可推得一种新的求积法，其步骤：(1) 取，由（6）式计算 再取，递推计算如下两步（）， （）对，计算 （3）对给定的误差，若有 ，即可停机，取，否则k+1代k转向（）。上述方法是逐次平分区间，是用复合梯形公式，在配合Richardson外推法，这种方法称为Romberg求积法。实际计算时如表：， 0123误差注（1）是复合梯形序列，是复合Simpson序列，是复合Cotes序列，他们与，统称为Romberg序列。（2）由于当较大时，再做组合意义不大，通常算到即可。（3），例 用Romberg 求积法计算，使误差不超过。解：将计算结果列表如下k010

39、.939793323 ，例 已知，求。 解：，4.4 Gauss型求积法前面提出的新的求积公式目的主要是改善截断误差。这节将讨论是否可以提高代数精度？在求积区间上节点数不变情况下是否可以通过选择节点，以提高求积公式的代数精度？一、Gauss型求积公式 1. 定义 设节点，求积公式， （1）其中 ，若（1）式至少有次代数精度，则称它为Gauss型求积公式，称为Gauss点，称为Gauss求积系数。2. Gauss点的性质定理 是Gauss点的充要条件：对次数不超过的任意多项式，与在上正交其中。证明：是Gauss点，任取次数不超过的多项式，是不高于的多项式，于是。若上式对任意次数不超过的多项式成立

40、，考虑次数不超过的多项式，记，次数不超过，=，而 ，即说明是Gauss点。3. Gauss-Legendre求积公式，是-1,1上以为权的正交多项式系。 设在-1,1上零点为，则是求积公式 （4）的Gauss点。例 （1）取的零点为节点，作求积公式令它对严格成立，得，从而得零阶Gauss-Legendre求积公式 （2）取的零点为节点作求积公式令它对精确成立：于是一阶Gauss-Legendre求积公式它有3次代数精度。二 Gauss型求积公式的截断误差、代数精度、数值稳定性和收敛性1. -,；2 . 代数精度为；Gauss型求积公式是数值稳定的 证明：先证，再证 设的实际计算值为，记，。证毕

41、 4. Gauss型求积公式是收敛的。例1 证明中矩形求积公式是高斯型求积公式 。 设充分光滑，试确定其余项。证明：因为 当时，左边=右边； 当时，左边=右边； 当时，左边=右边，所以上求积公式有1次代数精度，是高斯型求积公式。=。 例2 确定，使下面求积公式的代数精度尽可能高并求其代数精度 。解：当时，左边=，右边=； 当时，左边=0，右边=；当时，左边=，右边=令 上面左边=右边得当时，左边=右边；当时，左边=右边=，所以代数精度为3次。 例3 确定，使下面求积公式的代数精度尽可能高并求其代数精度 解：当时，左边=右边； 当时，左边=0，右边=；当时，左边=，右边=；令 上面左边=右边得第

42、五章 常微分方程数值解法 5.1 引 言回顾：微分方程的分类？能求解哪些微分方程？是否所有微分方程都能求解？不能求解的微分方程怎么办？ 以如下初值问题为例来介绍微分方程的数值解法。讨论初值问题 （1）有解的条件：设 f (x, y)在区域 R ：a x b, y 上连续，对y满足Lipschitz条件：其中，L 为正常数。设精确解为。所谓数值解法就是求在一列点a = x0 x1 x2 x n = b上的近似值，。通常取 xi = a + ih , h=(b a) / n 。 Euler法和改进的Euler法 Euler 法1、Euler公式 因 ，对 所以 ，将 y (xi)用差商近似： ，

43、再用yi近似代替y (xi)，由此得Euler公式(h = xi+1 xi): （2）注：Euler法的几何意义是 用过点的一条折线近似代替过点的积分曲线。例1 用Euler法求初值问题: 解 （1） Euler公式：，这里，=（2） 算得 y1 = y(0.1) = 0, y2 = y(0.2) = 0.01, y3 = y(0.3) = 0.029, 2、Euler法的截断误差 局部截断误差和整体截断误差（1）局部截断误差假设在计算时用到的是精确值即时产生的误差记 ，将在处泰勒展开 （3）称为局部截断误差。（2）整体截断误差= =+ 反复利用可得3、Euler方法的数值稳定性 仅就 （4）

44、讨论之。 记是（的）实际计算解， 假设用某个数值方法求解常微分方程，对给定，产生误差，若这个误差在计算后面的，中所引起的误差按绝对值均不增大，则称这个数值方法对和复常数是绝对稳定的。绝对稳定区域。考虑Euler方法 若 ，误差就不增大，这时Euler方法绝对稳定。Euler方法的绝对稳定域为：。 改进的Euler法 因 ，对 所以 ，将用差商近似：再用yi近似代替y (xi)，由此得向后的Euler公式(h = xi+1 xi): （5）其几何意义为过(xi , yi)作射线与过(xi+1 , y(xi+1)点处曲线的切线平行。结合（2）和（5）式，得梯形公式: （6）例2 用梯形法解例1。解

45、 公式：yi+1 = yi + 0.05 (xi yi + xi+1 yi+1), 或 yi+1yixi + 0.005)/1.05, 解得 y1 = y(0.1) = (0+0+0.005)/1.05 = 0.004762, y2 = y0.1+0.005)/1.05 = 0.01859, 余见下表.xiEuler法向后Euler法改进的Euler法预测-校正法k1k20000000910Euler预测校正法梯形公式是隐式的，求yi+1必须解一个方程。一般我们采用迭代法求解：用Euler公式求出初始近似 yi+1(0)，然后用梯形公式进行校正： 直到满足 | yi+1(k+1) yi+1(k

46、) | ， 取yi+1= yi+1(k+1) ，再转到下一步计算。只要f (x, y)对y满足Lipschitz条件，当h足够小时迭代序列收敛。实际应用梯形公式公式时，只需进行一次迭代就可以了，由此得梯形公式法的另一形式，或叫Euler预测校正法。 （7）写为下面的形式更便于上机计算： （8）例3 用Euler预测-校正法解例1(y = x y, y(0) = 0)。解 x1 ：k1=00=0，k2，y1；x2：k1，k2， y2(0.095+； 5.3 RungeKutta 法 RungeKutta法1、由微分中值定理 ， （1）Euler公式与上式比较：Euler法是用点的斜率近似代替得出

47、的。上式与Euler预测校正法比较：Euler预测校正法是用、上的斜率的平均值近似代替得出的。启示：当取斜率的点数确定时，可对点的位置和平均权引进待定参数，然后选择参数已取得较高的局部截断误差的阶数。这就是RungeKutta方法的精髓所在。考虑： （2）选择参数，使 （3）的阶尽可能高，其中设有适当的光滑性。因为 y (x) = f (x, y (x), y (x) = f x + f y y (x) = f x + f y f , 所以 （4） （5）将上两式代入（3）式得 （6）令 （7）代入（6）式 如取，得，代入（2）式即得预估-校正公式； 如取，得，代入（2）式即得 （8）称之为中

48、点公式。二级二阶R-K公式。 是否可选择使（对）的阶更高？可以证明这是不可能的。 类似地可以推导三级、四级以至一般的m-级R-K公式 ， ， ， ， 可推两个主要公式1、四级四阶经典R-K公式 (8)2、四级四阶Gill公式四阶RongeKutta公式不太复杂，精确度也够高。高于四阶的公式计算量大，精确度也并不一定提高，有时还会降低。 RK法的绝对稳定区域以四级四阶经典RungeKutta方法为例针对讨论，得=于是 =相应的误差方程=由此得绝对稳定区域5.4 单步法的收敛性计算公式访法的收敛性：对，有 定理 若（1）（2）在上满足，（3）则 （1）；（2）当时，方法是收敛的。5.5 解椭圆型方

49、程的差分法设有Poisson方程第一边值问题：Laplace算子 区域的离散化， ，步长；节点；相邻节点；内部节点；边界节点； 微分方程的离散化 （3） （4）引入网络差分算子+ （5） （6）称为微分方程（1）的五点差分格式，他的解称为差分解。注 （1）正方形网格，（6）式变为 （7） （2）用（6）逼近（1）的截断误差 有 称差分方程（6）与微分方程（1）是相容的。 边界条件的处理1、矩形区域（矩形网络），2、一般区域 （1）直接转移：若边界点S正好在上，则 若边界点P不在上，可在上取与P最靠近的网格上的点R，用代： （2）线性插值 边界点P也可取B、C两点或R、Q两点的线性插值作为。 设

50、，取或边界点的方程与内点的方程构成一个方程组。其方程个数与未知量的个数相等常用矩形域D：A是不可约对角占优，且是大型稀疏带状矩阵。采用迭代法求解。例2 用标准4阶Ronge-Kutta法解初值问题 解 这里, x0 = 0, y0 = 1, h = 0.2. 利用标准4阶Ronge-Kutta公式有计算结果如下表. ixiyik1k2k3k400111234517 (3) 常微分方程组和二阶常微分方程初值问题的RongeKutta公式设有微分方程组初值问题 (9)其求解公式为: (10)以上公式写为向量形式与单个方程的公式形式上完全相同, 证明也相仿. 其几何意义可理解为参数方程表示的曲线.对

51、二阶方程 (11)设 x= y, 则 .y = f (t, x, y). 在(10)式中取 g (t, x, y) y, 得 (12)5.6 Adams 法 RongeKutta法优点: 是单步法, 由前一步的 yi 值就可以推得后一步的yi+1 值. 缺点: 当 f 较复杂时计算量较大. 未利用前面几步算得的 y 值.我们知道, 初值问题 y = f (x, y), y (x0) = y0等价于积分方程 过前几步算得的(xi , yi)作 f (x, y (x)的插值多项式pn(x), 设余项为R n(x), 则 f (x, y) = pn(x) + R n(x) 舍去余项后得到线性多步法公

52、式: (1) 5.6.1 Adams 显式取q +1个结点 xi , xi1 , , xiq ,并作Newton后差插值多项式, 则 (2)其中 x = xi + sh , fi = f (xi , yi). 将(2)式代入(1)式, 得 (3)这里, . (3)式即为Adams显式.由此得Rq = h q+2 y(q+2) () q . (x iq, x i) (4) m是多项式积分, 列表如后. 通常将(3)中的 m f 用各点已知函数表示, 如 q=2, 3 时有公式 (5) (6)Adams法系数 m 、m 表m01234 m11/25/123/8251/720m11/21/121/2

53、419/7205.6.2 Adams隐式 显式公式作插值多项式时未包括xi+1, 精确度较差. 若取q +1个结点 xi+1 , xi , , xiq+1 ,并作Newton后差插值多项式, 则 (7)其中 x = xi+1 + sh .将(7)式代入(1)式, 得 (8)这里, , m = 0, 1, 2, q. (8)式即为Adams隐式.类似于Adams显式的余项求法, 可得Adams隐式的余项为: Rq = h q+2 y(q+2) q+1, (xiq+1, xi+1) (9) m 的计算结果如节表.将 mfi 用各点已知函数表示, q = 2, 3 时有公式: (10) (11)5.

54、6.3 Adams预测校正公式通常将Adams隐式和显式联立使用, 可得Adams预测校正公式. 以q = 3为例: (12)用pi与ci依次表示第i步yi的预测值与校正值, 精确值仍记为yi , 由误差估计式(4), (9)可得 相减得 从而有 由此得 由此给出预测校正法(12)的另外一种模式(简称为PMCM模式):(13)开始时可取 pi = ci = 0 , 以后每步按上式计算. 因为Adams法是多步法, 不能自动开始, 前几个值依赖于其它方法, 例如我们可以:(1)用RongeKutta法求出开始值.(2)使用 y (x) 的Tailor展开式.5.7 二阶线性常微分方程边值问题在解

55、数学物理问题时经常遇到各种边值问题. 这里我们只介绍二阶线性常微分方程边值问题 (1)其中 p (x), q (x), f (x)为区间a, b上足够光滑的已知函数且q (x) 0, , 为已知常数.将区间 a, b 分成n等分, 将各结点导数用差商表示, 可化为差分方程组. 设y , y 的误差均为O (h 2). 设 pi = p (xi), qi = q (xi), fi = f (x i), 代入方程(1)得: 两边乘以h 2, 整理得 (2)即 (3)例 用差分法解边值问题 解 这里，h = 0.1, x i i, p = 0, q = 1, f (x) = x, i = 0, 1,

56、 , 10. 代入(3)式得ixi差分解yi精确解yi123456789解得结果见右表.最后一列为精确解 2 sh x / sh1x.关于差分方程的收敛性有如下的收敛定理 设 2 pi h 2, qi 0, i = 1, 2, , n 1, 则二阶差分方程(2)有唯一解, 且当h 0时收敛于微分方程(1)的准确解.第六章 非线性方程求根 掌握非线性方程求根的每一种方法的基本思想、理论依据和基本步骤6.0 引 言1、问题 设，则至少存在，使得，如何求？我们都会解一元一次方程和一元二次方程，三次以上方程就不会解了。事实上，三次和四次方程的求根公式很复杂，五次以上代数方程一般无求根公式，超越方程也只

57、能求近似解。2、求近似解的步骤(1) 判定根的存在性，根据连续函数介值定理 f C a, b, f (a) f (b) 0 (a, b)，使得 f ( ) = 0(2) 确定根的分布区间: 可列函数表或求导数确定单调区间。(3) 根的精确化缩小区间。3、常见的求近似解的方法：二分法、牛顿法、弦割法和抛物线法等 6.1 二分法基本思想：对有根区间不断进行二分；理论依据：区间套定理和收敛序列的性质； 二分法的思想：逐渐二分有根区间，得一系列有根区间 ，用某个区间的中点 作为根的近似值。二分法的理论依据：得三个序列 ，有，且 二分法的算法（自己完成）：输入，计算；判别，是输出，否判别，是，否判别，是

58、输出，停机；否转（2）结束例1 用二分法求方程 x 2 2 = 0的根。解 设f (x) = x 2 2，因f (1) = 1 0, 故在（1, 2）内有根。取 x1， 得 f (1.5) = 1.25 0。 取 x2，得 f (1.25) = 0.43 0。令 ，利用介值定理可证至少有一根，再利用反证法证只有唯一根（即）。2) ，由条件（2）知，设，同理可得，由数学归纳法得对，由（3）式得xk。由Lipschitz条件 ；将 代入上式得。 如何构造等价方程？（不同的构造会得到不同的收敛性） 把方程 f (x) = 0改写为x = g (x)可以有多种方案。例1 x3 x 1 = 0 在附近有

59、根，把它写成不同等价形式，判断各迭代公式的收敛性并迭代求解(保留6位有效数字，下文同，不再声明)。 x = x3 1； 。解 设g1 (x) = x3 1，g1 (x) = 3 x2，| g1 (1.3) | = 3.4 1。故迭代序列xk+1 = xk3 1发散。 设。故迭代序列收敛。用它进行迭代，得结果如下表。因此，x 。k012345678xk例 用简单迭代法求在2，3上的根。为求在2，3上的根，可构造三个同解方程（1），（2），（3）得三个迭代： 一个收敛快；一个收敛慢；一个发散。 （1），（2），（3）取，结果如下表0222112-53-1254-19530055 迭代法的局部收敛性

60、若 ，使 对都收敛，则称上迭代是局部收敛的。定理2 若在根的某个邻域内有连续导数，且，则上迭代有局部收敛性。 收敛阶定义 设某迭代有局部收敛性，记ek = s xk 0, k = 0, 1, , 如果存在正实数r和c使成立，则称该迭代是r阶收敛的。r=1时又称有线性收敛速度；r=2称平方收敛。例 xk k，ekk 有线性收敛速度(每算一步小数点后有效数字每次增加1位)。 yk 3 k 有三阶收敛速度(小数点后有效数字成三倍增加)。定理3 设I = a, b， ，s = g (s) I，g (k) (s) = 0, k = 1, 2, m 1, g (m) (s) 0，则对任何x0 I (x0