基于matlab的数值计算中的优化技术

1、Hefei University 毕业论文（设计） BACH ELOR DI SSERTATI ON 论文题目： 基于 matlab 的数值计算中的优化技术 学位类别： 理 学 学 位 学科专业： 信 息 与 计 算 科 学 基于 matlab 的数值计算中的优化技术 中文摘要 优化是人们寻求的目标，数值计算中优化技术采用的好，能从时间与空间上得到 巨大的好处。一个算法除了正确外，还要空间能存贮程序数据，且运行时间短，因而 算法优化技术就很重要，一个程序无法存贮到计算机内存中，或运行时慢得无法等候 是没有任何实际意义的。由于数值计算的优化技术有很多方面，为此选用数值积分进 行说明，数值积分计算

2、有很多的方法，用不同的方法所计算的积分精确度不同，所需 要的时间也不同，通过一些实例的分析，对优化技术进行归纳与总结。 矚慫润厲钐瘗睞 枥庑 赖。 本论文是基于 matlab 在数值计算中的优化技术，优化技术是算法设计的重要而 关键的课题，本论文选取数值分析中的一些著名的优化技术进行讨论，并在 matlab 中加以实现，通过 tic 、tuc 、 cputime 等函数的使用对其进入深入分析。 聞創沟燴鐺險 爱 氇谴净。 关键字：数值计算；优化技术； matlab ；数值积分 Numerical optimization technique based on MATLAB ABSTRACT O

3、ptimization is that people seek target and numerical optimization technique used is good. It can get huge benefits from the time and space. An algorithmnot only in 残骛 order to right, but also space can store data, and short running time. So the algorithm optimization technique is very important and

4、a program cannot be stored in the computer memory, or run slow and practical of significance. The optimization technology of numerical calculation is used in many respects. The numerical integration is described and numerical integral calculation has many methods. Integral accuracy is calculated by

5、different methods and the time required is different also, and through the analysis of some examples, and summarizes the optimization technique. 楼諍锩瀨濟溆塹籟。 This paper is to optimize the technology based on MATLAB in the numerical calculation. Optimization technique is the key task in the algorithm de

6、sign, this paper selects some well-known optimization techniques in numerical analysis are discussed, and implemented in MATLAB, through in-depth analysis using tic, tuc and cputime and other functions of the entry. 酽锕极 額閉镇桧 猪訣 锥。 KEYWORD: numerical calculation; optimization techniques;MATLAB;numeri

7、cal integration 彈贸摄 尔霁毙 攬砖卤庑。 第一章 前言 1 謀荞抟箧飆鐸怼类蒋薔 第二章 数值积分的计算 2 厦礴恳蹒骈時盡继價骚 2.1 数值求积公式的构造 2 茕桢广 鳓鯡选块 网羈泪 2.1.1 求积公式的推导 2 鹅娅尽 損鹌惨歷 茏鴛賴 2.1.2 几个低次牛顿 - 科特斯求积公式 4 籟丛妈 羥为贍偾 蛏练淨 2.2 复化求积公式 6 預頌圣 鉉儐歲龈 讶骅籴 2.2.1 复化梯形求积公式 6 渗釤呛 俨匀谔鱉 调硯錦 2.2.2 复化辛浦生求积公式 7 铙誅卧 泻噦圣骋 贶頂廡 2.2.3 复化科特斯求积公式 7 擁締凤 袜备訊顎 轮烂蔷 2.3 高精度数值积分

8、算法 8 贓熱俣 阃歲匱阊 邺镓騷 2.3.1 龙贝格求积公式 8 坛摶乡 囂忏蒌鍥 铃氈淚 第三章 线性方程组的求解 10 蜡變黲癟報伥铉锚鈰赘 3.1 线性方程组的介绍 3.2 线性方程组的迭代法 10 買鲷鴯 譖昙膚遙 闫撷凄 11 綾镝鯛 駕櫬鹕踪 韦辚糴 3.2.1 Jacobi 迭代法 11 驅踬髏 彦浃绥譎 饴憂锦 3.2.2 Gauss-Seidel迭代法 12 猫虿驢 绘燈鮒诛 髅貺庑 3.2.3 SOR 迭代法 13 锹籁饗 迳琐筆襖 鸥娅薔 14 構 氽頑黉碩饨 荠龈话骛 第四章 各种求积公式的 MATLAB编程实现与应用 4.1 对数值积分运行结果及其分析 14 輒峄陽

9、 檉簖疖網 儂號泶 4.1.1数值积分运行结果 14 尧侧閆 繭絳闕绚 勵蜆贅 4.1.1数值积分运行结果分析 15 识饒鎂 錕缢灩筧 嚌俨淒 4.1 线性方程组运行结果及其分析 15 凍鈹鋨 劳臘锴痫 婦胫籴 4.2.1线性方程组运行结果 15 恥諤銪 灭萦欢煬 鞏鹜錦 4.2.2线性方程组运行结果分析 16 鯊腎鑰 诎褳鉀沩 懼統庫 附录 17 硕癘鄴 颃诌攆檸 攜驤蔹 参考文献 23 阌擻輳嬪諫迁择楨秘騖 致谢 24 氬 嚕躑竄贸恳 彈瀘颔澩 第一章 前言 数值计算是有效使用数字计算机求数学问题近似解的方法与过程，以及由相关理 论构成的学科。数值计算主要研究如何利用计算机更好的解决各种数

10、学问题，包括连 续系统离散化和离散形方程的求解，并考虑误差、收敛性和稳定性等问题。从数学类 型分，数值运算的研究领域包括数值逼近、数值微分和数值积分、数值代数、最优化 方法、常微分方程数值解法、积分方程数值解法、偏微分方程数值解法、计算几何、 计算概率统计等。随着计算机的广泛应用和发展，许多计算领域的问题，如计算物理、 计算力学、计算化学、计算经济学等都可归结为数值计算问题。 釷鹆資贏車贖孙滅獅赘。 优化是人们寻求的目标，数值计算中优化技术采用的好，能从时间与空间上得到 巨大的好处。一个算法除了正确外，还要空间能存贮程序数据，且运行时间短，因而 算法 优化技术就很重要， 一个程序无法存贮到计算

11、机内存中， 或运行时慢得无法等候 是没有任何实际意义的。 在研究基于 matlab 在数值计算中的优化技术有很多方面求 数值积分就是具有代表性的一点。 怂 阐譜鯪迳導嘯畫長凉。 求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出 来，因此能够借助微积分学的牛顿 - 莱布尼兹公式计算定积分的机会是不多的。另外， 许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的 定积分，也不能用不定积分方法求解。由于以上原因，数值积分的理论与方法一直是 计算数学研究的基本课题。通过这个课题的研究，我们将会更好地掌握运用数值积分 算法求特殊积分函数的定积分的一些基本方法

12、、理论基础；并且通过matlab 软件编程 的实现，得出计算数值积分的最优化的方法，并应用于实际生活中。谚辞調担鈧谄动禪泻 類 第二章 数值积分的计算 2.1 数值求积公式的构造 人们根据积分的定义得出 Newton-Leibniz 求定积分的公式 , 但是这些公式并不是 能求出所有式子的积分 , 而是针对许多特殊的例子 , 但是有许多都是球不出来的 , 如 : 1 sin x 2sin x esin x sin xdx,dx 等。 嘰觐詿缧铴嗫偽純铪锩 。 01 所以采用积分的几何意义来设计出积分公式从而求出近似值。 2.1.1 求积公式的推导 建立数值积分公式的途径比较多，其中最常用的优两

13、种： 1)对于连续函数 (x) ，优积分中值定理： f (x)dx (b a) f ( ) ( a,b) (2.1.1) 其中 f ( )是被积函数 f (x) 在积分区间上的平均值。因此，如果能给出求平均值 f ( ) 的一种近似方法，相应地就可以得到计算定积分的一种数值方法。 熒绐譏钲鏌觶鷹緇機 库。 (2)先用某种简单函数 (x)近似逼近 f(x) ，然后 (x)在 a,b区间的积分值近似 表示 f (x)在 a,b区间上的定积分，即取 bb f(x)dx (x)dx (2.1.2) aa 一般情况下，我 们可以 取 (x)为前面介 绍的 f(x) 插值多项式 L(x)或拟合多项式 P(

14、x) 进行近似计算。若取( x)为插值多项式 Ln ( x) ，则相应得到的数字微分公式就是插 值型求积公式 鶼渍螻偉 阅劍鲰腎邏蘞。 把区间 a,b n等分，其分点为 xi a ih(i 0,1,2,L ,n) 、h ba ，过这 n 1个节点， 可以构造一个 n 次插值多项式： Ln(x) w(x) i 0 (x xi )w (xi) f (xi ) (2.1.3) 其中 w(x) (x x0)(x x1)L (x xn)，用 Ln( x)代替被积分函数 f (x) 则有 f (x)dx b a Ln (x)dx w(x) i 0 (x xi)w (xi) f (xi ) dx 其中 Ai

15、 b a(x w(x)dx xi )w(xi) n i0 b a (x Aif (xi) i0 上式叫做牛顿 - 科茨公式，使用牛顿 x a th ，于是 w(x)dx xi)w(xi) f (xi ) (2.1.4) - 科茨公式关键是计算系数 Ai ，用变量替换 n1 w(x) w(a th) hn 1t(t 1)L (t n) (2.1.5) 而w(xi) hn( 1)n i (i!)(n i)! 这样 Ai b a (x w(x)dx xi)w (xi) i)hdx nhn 1t(t 1)L (t n) 0 ( 1)n i hn(i!)(n i)! h(t (2.1.6) ( 1)ni

16、h nt(t 1)L (t n)dx i!(n i)! 0 (t i) (2.1.7) (2.1.8) C(n) ( 1)n ih n t(t 1)L (t n)dx Cii!(n i)! 0 (t i) dx Ai (b a)Ci(n) 这时 Ci 是不依赖于函数 f (x) 和区间 a,b 的常数，可以事先计算出来，叫做牛顿 科茨系数。 表 2.1 牛顿 - 科茨系数 n Ci(n) 1 1 1 2 2 2 1 4 1 6 6 6 3 1 3 3 1 8 8 8 8 4 7 16 2 16 7 90 45 15 45 90 5 19 25 25 25 25 19 288 96 144 14

17、4 96 288 6 41 9 9 34 9 9 41 840 35 280 105 280 35 840 7 751 3577 1323 2989 2989 1323 3577 751 17280 17280 17280 17280 17280 17280 17280 17280 8 989 5888 928 10496 4540 10496 5888 989 28350 28350 28350 28350 28350 28350 28350 28350 2.1.2 几个低次牛顿 - 科特斯求积公式 1、梯形求积公式 定义 2.1 在牛顿 - 科特斯求积公式中，如果取 n 1，用一次多项式代

18、替被积函数， 即用梯形面积代替曲边梯形的面积，则有 纣忧蔣氳頑莶驅藥悯骛。 f (x)dx L1( x)dx (b a)c0(1) f (x0) c1(1) f (x1) (1) (1) 1 其中， f(x0) f (a), f(x1) f(b) 查表可得 c0(1) c1(1)代入上式得出 颖刍莖蛺饽亿顿裊赔 2 泷。 f(x)dx L1(x)dx (b a) 2 f (a) f (b) (2.1.9) 称式 (2.1.9) 为梯形求积公式 根据牛顿 - 科特斯求积公式的误差理论，梯形求积公式的误差估计为 R1(f) (2) ( ) (2)! b a(x a)(x b)dx (b12a)3f

19、() 12 f ( )是被积函数 f (x)二阶导数在 点的取值， a,b 濫驂膽閉驟羥闈詔寢賻。 2、辛浦生求积公式 定义 2.2 在牛顿 - 科特斯求积公式中， 如果取 n 2 ，用二次多项式代替被积函数， 即曲边用抛物线代替，则有 銚銻縵哜鳗鸿锓謎諏 涼。 bb f (x)dxL2(x)dx (b aa a)c0(2) f (x0) (2) c1(2) f (x1) c2(2) f(x2) 其中， x0 ab a， x1 2 ， x2 b, 查表可得 c0(2) c(22) 1 ， c1(2) 2 ，代入上式得出 挤 63 貼綬电麥结鈺贖哓类。 b b 1 2 a f (x)dx a L

20、2(x)dx (b a)6 f(a) 3 f ( a2b) 1 16f(b) 2.1.10 称式 2.1.10 为辛浦生求积公式，也称抛物线求积公式。 辛浦生求积公式的误差估计式 R2(f ) 41! (b120a) 5 f (4) ( ) (b a)5 f (4)( ) 2880 f ( ) a,b 3、科特斯求积公式 定义 2.4 在牛顿 - 科特斯求积公式中，如果取 n 4 时，牛顿科特斯公式为 赔荊紳 谘侖驟辽輩袜錈。 f (x)dx ba 90 7f (x0) 32f (x1) 12f(x2) 32f(x3) 7f (x4) 2.1.11 称式 2.1.11 为科特斯求积公式 同理可

21、求得其误差估计式 R4 f 2(b a)(b a)6 f(6)( )( (a,b) 945 4 2.2 复化求积公式 2.2.1 复化梯形求积公式 在上一节求积分的过程只是求粗约的近似值，所以应根据积分的可加性，可以将 区间分为许多部分使得积分值更加接近精确值，从而优化了梯形积分公式， 辛普生积 分公式和科特斯积分公式，这就是复化求积分公式的思想。塤礙 籟馐决穩賽釙冊庫。 ba 定义 2.5 将积分区间 a,b 进行 N 等分，记为 h b a, 裊樣祕廬廂颤谚鍘羋蔺。 N xk a kh 在每个小区间 xk ,xk 1 (k 0,1,L N 1) 上用梯形公式求和，得 仓嫗盤 紲嘱 珑 詁鍬

22、齊驁。 f (x)dx N1 k0 xk 1 k1 f (x)dx xk N 1 h hf(xk) f(xk 1) k 0 2 若将所得的近似值记为 TN ，整理得 f (x)dx h h2f(a) N1 f(b) 2f (xk) TN 2.2.1 k1 称式 2.2.1 为复化梯形公式。记为TN 复化梯形公式的截断误差 RT(f) (b a)h2 f 12 2.2.2 复化辛浦生求积公式 在辛普生积分公式上加以复化可以得到复化辛普生积分公式 定义 2.6 将积分区间 a,b 分成 N 2m等分，分点为 xk a kh， (k 0,1,L ,N 1) h b a 在每个小区间 x2k 2,x2

23、k(k 1,2,L ,N 1)上。用 Simpson 公式求积分， 绽萬璉轆娛閬 N 蛏鬮绾瀧。 得到 SN f (x)dx a h h3f(a) f (b) m1 2 f (x2k ) k1 m 4 f (x2k 1) k1 2.2.2 式 2.2.2 就称为复化辛浦生求积公式。记为SN 如果 f(x) C (4) a, b ， 则由 Simpson 插值余项公式可得复化公式的截断误差为 骁顾 燁鶚巯瀆蕪領鲡赙。 b h m 1 m m (2h)5 RS(f)a f(x)dxhf(a)f(b) 2f(x2k)4f(x2k1)(2h) f(4)() a 3 k1 k1 k 1 2880 x2k

24、 2,x2k 2.2.3 复化科特斯求积公式 定义 2.7 将积分区间a,b等分为 N个子区间x4k 4,x4k，每个子区间的中点 x4k 2, (k 0,1,L ,N 1) ，子区间长度 h b a， 在每个子区间上用科特斯公式求和， N 得 瑣钋濺暧惲锟缟馭篩凉。 N 12 f(x4k 2) k1 hN f (x)dx h 7 f (a) 32 f (x4k 3) 90 k 1 N1 32 f (x4k 1) 14 f(x4k) 7f(b) 2.2.3 k1 k1 式 2.2.3 就称为复化科特斯求积公式，式中 h b a ，xk a kh (k 0,1,2,L ,2N 1) N4 类似地

25、可以推出复化科特斯公式的截断误差为 R4(N)(f)2(b a) ( h)6 f (6) ( ) ( (a,b) 945 4 2.3 高精度数值积分算法 求积分时，复化积分公式采用逐步分段法，是一种比较有效的方法，但是它也存 在许多的弊端，它收敛于积分真值的速度缓慢，从而人们在复化求积分公式上进行改 进。在求积分时步长的大小 也会影响到积分的效果，步长太长积分值就不会太精确， 步长太短则会增加许多的运算量。运算的时间也会随之增加。在计算器中编程计算结 果要花费大量的时间。以下采用变步长的计算方法，从而避免了这一点。 鎦诗涇艳 损楼紲 鯗餳類。 2.3.1 龙贝格求积公式 梯形法的算法简单，单精

26、度低，收敛的速度缓慢。由此引出了龙贝格公式。由梯 形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果T2N 作为积分 I 的近似值，其误差近似值为1(T2N TN ) 。可以设想，如果用这个误差作为T2N 的一种补 3 偿，即将 栉缏歐锄棗鈕种 鵑瑶锬。 13(T2N TN)= 4T24N 1TN 作为积分的近似值，可望提高其精确度。 直接根据复化求积公式，不难验证 SN T2N 13(T2N 4T2NTN 41 (2.3.1) 这说明，将区间对分前后两次复化梯形公式的值，作线性组合恰好等于复化辛浦 生公式的值 SN，它比 T2 N更接近于近似值。 同样，用 S2N于SN作线性组合

27、会得到比 S2N更 精确的值，且通过直接验证可得 辔烨棟剛殓攬 瑤丽阄应。 1 CNS2N(S2NSN ) 15 4 S2NSN 42 1 (2.3.2) 用C2N与 CN作线性组合，又可得到比C2N更精确的值，通常记为 RN ,即 (2.3.3) C 1 (C C ) 4 C2NCN C2N(C2NCN )3 63 4 1 式 (2.3.3) 就称为龙贝格求积公式。 上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推方法。 我们将序列 TN , SN CN 和 RN 分别称为梯形序列、辛浦生序列、科特斯序列和龙 贝格序列。由龙贝格序列当然还可以继续进行外推，得到新的求积序列。峴

28、扬斕滾 澗辐滠 兴渙藺。 在积分区间逐次分半的过程中，利用外推法算式将粗糙近似值 Tn 逐步加工成越来 越精确的近似值 Sn,Cn,Rn L 。也就是说，将收敛速度缓慢的梯形序列 T 2k 逐步加工成 收敛速度越来越快的序列 S2k, C2k , R2kL 。根据这种原理设计的计算积分近似值的方 法称为龙贝格积分方法，又称为数值积分逐次分半加速收敛法。 詩叁撻訥烬忧毀厉鋨骜。 利用龙贝格序列求积的算法称为龙贝格算法。这种算法具有占用内存少、精确度 高的优点。因此，成为实际中常用的求积方法。在优化技术方面的考虑龙贝格方法比 则鯤愜韋瘓賈晖园栋泷。 较合适的选择，在第四章会用程序进行讨论说明 第三

29、章 线性方程组的求解 上一章讲述了有关于优化技术在数值积分计算方法的应用，为了更加体现优化技 术的应用本章将会讨论优化技术在线性方程组中的应用。实际中，存在大量的解线性 方程组的问题。很多数值方法到最后也会涉及到线性方程组的求解问题：如样条插值 的 M和 m 关系式，曲线拟合的法方程，方程组的Newton 迭代等问题。 胀鏝彈奥秘孫戶 孪钇 賻。 求解线性方程组有很多的方法， 如 gauss 消去法， 按比例主元消去法， 用 Cholesky 分解解线性方程组，平方根法和追赶法等等。 鳃躋峽祷紉诵帮废掃減 。 3.1 线性方程组的介绍 般地设 n 阶线性方程组为 表示成矩阵形式 其中 a11

30、a12 Aaij n n a21 a22 an1 an2 a11x1 a12x2 a1nxn a21x1 a22x2 a2nxn an1x1 an2x2 annxn Ax b x1 b1 ，x x2 ， b b2 xn bn bn ann a2n a1n b1 b2 A 为系数矩阵 高斯消元法是按照消元和回代两个过程。高斯消元法的改进为高斯主元消元法， 并且主元消元法主要有列主元，按比例主元和全主元。 稟虛嬪赈 维哜妝扩踴粜。 高斯消元法的基本思想： 首先将 A 化为上三角阵，再回代求解 10 0 an(22) an(2n) bn(2) 00 (3) an3 (3) ann bn(3) 类似下

31、去我们有 第k步 ai(kk) 第k行a(kik) akk 第i行,i k 1, ,n 1) 2) 3) Mn) (b1 (2b2 (3b3 M(b (1)1n (2)n2 (3)3n M(n aaa a LLLO ) 1 (a (2)a ) 1 (a b1 b2 b nn a1 a2 a 22 12 aa a 11 a1 a2 a 1 2 3 b1(b2 ( n 2n 3 a1(a2( 32 1( aa 12 (2 aa a11 0 a1 2 (b n1 步以后，我们可以得到变换后的矩阵为： (2 n2n 3 1(2 ( aa 2 (a 22 a1(a a1 0 an(nn) 从此再回代可以

32、 解出线性方程组的解。 但是高斯解线性方程组一般都是针对中小型的，一下介绍几种线性方程组的迭代法，从而求 解线性方程组的近似解，利用优化技术判断哪个方法最优。 陽簍埡鲑罷規呜旧岿錟。 3.2 线性方程组的迭代法 3.2.1 Jacobi 迭代法 设 n 阶线性方程组的系数A (aij )n n 非奇异( nonsigular )，且 aij 0(i 1,2,L ,n) 。将方程组改为 11 x1 a1 ( a11 a12x2 a13x3 L a1nxn b1), x2 a1 ( a21x1 a23x3 L a2nxn b2), a22 (3.2.1) L, 1 xn ( an1x1 an2x2

33、 L an,n 1xn bn). ann 任取 x(0) (x1(0) ,x2(0) ,L , xn(0) )T ，将各分量代入上式的右边得 (1) x1 a1 ( a11 (0) a12x2 (0) a13x3 L (0) a1nxn b1), (1) x2 a1 ( (0) a21x1 (0) a23x3 L (0) a2nxn b2), a22 L, (1) xn 1 ( (0) an1x1 an2x(20) L ax(0) an,n 1xn 1 bn) ann (3.2.2) 将 x(1) (x1(1), x(21) ,L , xn(1) )T代入上式的右边， (2) x1 a1 (

34、a11 (1) a12x2 a13x3 L (1) a1nxn b1), (2) x2 a1 ( (1) a21x1 (1) a23x3L (1) a2nxn b2), a22 L, xn(2) n 1 ( (1) an1x1 an2x2(1) L ax(1) an,n 1xn 1 bn). ann (3.2.3) 以此类推，可得 (k 1) a1 ( ann (k) an1x1 an2x2(k) ax(k) an,n 1xn 1 bn). 由此可得向量序列 x(k)(k 1,2,L ) 。称由迭代式建立的迭代法为 Jacobi 迭代法。 x1(k 1) 1 ( a11 (k) a12x2 (

35、k) a13x3 L a1nxn(k) b1), (k 1) x2 a1 ( (k) a21x1 a23x3(k) L a2nxn(k) b2), a22 L, 2 (3.2.4) 3.2.2 Gauss-Seidel 迭代法 在 Jacobi 迭代法中，每次迭代计算 x(k 1)时用的是前一次迭代的全部分量x(jk 1)(j 1,2,L ,n)。实 12 际上，在计算分量 x(jk 1)时，最新的分量 x1(k 1) ， x(2k 1)，L ， x(jk11)已经算出，但没有被利用，而且， 如果 Jacobi 迭代收敛，最新算出的分量xi(k 1)(j 1,2,L , j 1)一般比 xi(

36、k)(j 1,2,L , j 1)的精度 更高。因此，可以对 Jacobi 迭代法加以改进，即在迭代过程中，每个分量计算出来之后，计算一 下分量时就利用最新计算出的近似结果， 具体地， 即用新分量 xi(k 1)( j 1,2,L , j 1) 去替换 (3.2.4) 右端的各项中 xi(k)( j 1,2,L , j 1) ，可得新的迭代公式 (k 1) 1 (k) x1 ( a11 a12 x2 x(2k 1) 1( a22 (k) a21x1 (k 1) a1 ( ann (k) xn an1x1 (k) a13x3 L a1n (k) xn b1), (k) a23x3 L a2n (

37、k) xn(k) b2), 2 (3.2.5) L, an2x(2k) L an,n (k ) 1xn 1 bn). 沩氣嘮戇苌鑿鑿槠谔應。 此式称为 Gauss-Seidel 迭代法。 为了提高收敛速度，对 Gauss-Seidel 第 k 1 个近似解 钡嵐縣緱虜荣产涛團蔺。 迭代法进一步用 Gauss-Seidel 迭代公式计算得到 (k 1) 1 xi(bi aii i1 (k 1) aij xj j1 n (k) aij xj ) ,(i 1,2,L ,n) ji1 将前一步迭代值 xi(k) 与 Gauss-Seidel (k 1) 迭代值 xi 做加权平均，即 3.2.3 SOR

38、 迭代法 xi(k 1) (1 w)xi(k ) wxi ， i 1,2,L ,n 其中 w 是参数，整理得 xi ( k 1) (1 w)xi(k) i1 n (k 1) (k) aij xj ) aij xj j1 ji1 w (bi aii i 1,2,L ,n 此式称为松弛迭代法，其中参数w 为松弛因子 当 w 1 时，式称为超松弛法； 当 w 1 时，式称为低松弛法； 当 w 1 时，式就是 Gauss-Seidel 。 一般称这些方法为 SOR 方法。 13 第四章 各种求积公式的 MATLAB编程实现与应用 4.1 对数值积分运行结果及其分析 31 不同的方法计算出来的积分求值的

39、结果不同，所以一下针对 I 1dx 进行求解，运用不同 的方法得到近似值，余项以及运行程序所需要的时间，在精度不同的情况下，再进行分析判断出 哪个方法最优。 懨俠劑鈍触乐鹇烬觶騮。 4.1.1 数值积分运行结果 表 4.1.1 精度为 0.5 10 3 的积分公式运行结果 方法 n 近似值 时间 梯形求积公式 1 1.468592464645272 0.029394761280958 辛普生求积公式 1 1.468592464645272 0.020265210920251 科特斯求积公式 1 1.468592464645272 0.012255210320241 复化梯形求积公式 8 1.4

40、64104033541694 0.141131524566874 复化辛普生求积公式 2 1.46410919941527 0.12400000000000 复化科特斯求积公式 2 1.46410919941527 0.12400000000000 龙贝格积分公式 2 1.464106276346608 0.002670 表 4.1.2 精度为 0.5 10 4 的积分公式运行结果 方法 n 近似值 时间 梯形求积公式 1 1.468592464645272 0.029394761280958 辛普生求积公式 1 1.468592464645272 0.020265210920251 科特斯求

41、积公式 1 1.468592464645272 0.012255210320241 复化梯形求积公式 8 1.464104033541694 0.148061234217113 复化辛普生求积公式 8 1.464104033541694 0.143152625817205 复化科特斯求积公式 4 1.46410180228421 0.28100000000000 龙贝格积分公式 4 表 4.1.3 精度为 0.5 10 5 的积分公式运行结果 方法 n 近似值 时间 梯形求积公式 1 1.468592464645272 0.029394761280958 辛普生求积公式 1 1.4685924

42、64645272 0.020265210920251 科特斯求积公式 1 1.468592464645272 0.012255210320241 复化梯形求积公式 16 1.464101769539943 0.277173394958087 复化辛普生求积公式 16 1.464101769539943 0.283381996592803 复化科特斯求积公式 8 1.46410161860649 0.56100000000000 龙贝格积分公式 14 表 4.1.4 精度为 0.5 10 6 的积分公式运行结果 方法 n 近似值 时间 梯形求积公式 1 1.468592464645272 0.0

43、29394761280958 辛普生求积公式 1 1.468592464645272 0.020265210920251 科特斯求积公式 1 1.468592464645272 0.012255210320241 复化梯形求积公式 32 1.464101624841357 0.567891012749313 复化辛普生求积公式 32 1.464101624841357 0.563100795217624 复化科特斯求积公式 8 1.46410161860649 0.59300000000000 龙贝格积分公式 表 4.1.5 精度为 0.5 10 7 的积分公式运行结果 方法 n 近似值 时间

44、 梯形求积公式 1 1.468592464645272 0.029394761280958 辛普生求积公式 1 1.468592464645272 0.020265210920251 科特斯求积公式 1 1.468592464645272 0.012255210320241 复化梯形求积公式 64 1.464101615745077 1.129131380037762 复化辛普生求积公式 32 1.464101615745077 1.153032881910490 复化科特斯求积公式 16 1.46410161519479 1.17000000000000 龙贝格积分公式 4.1.1 数值积分

45、运行结果分析 根据 matlab 积分的命令的准确值为 I=1.464101615137754 。从运行的结果可以看 出梯形积分公式和辛普生积分公式运行的结果都差不多，但是科特斯公式运行的结果 较为准确，但是运行的时间相比科特斯积分公式的时间比较长。相同的精度但是科特 斯所分割的节点少。 謾饱兗争詣繚鮐癞别瀘。 这是因为梯形公式、 辛普生公式是低精度公式， 但对被积函数的光滑性要求不高， 他对对被积分光滑性较差的积分很有效。特别是梯形积分公式对被积分函数式周期函 数积分时，效果更加突出。高阶科特斯求积分公式稳定性差，收敛较慢。从而，为了 提高收敛速度建立的复化梯形积分公式，复化辛普生积分公式。

46、但是相比之下，龙贝 格积分公式是算法简单，是一个很好的加速方法。从优化技术的角度看，一般都是选 择龙贝格积分方法。 呙铉們欤谦鸪饺竞荡赚。 4.1 线性方程组运行结果及其分析 4.2.1 线性方程组运行结果 线性方程组 15 10 x1 x2 2x3 7.2 x1 10 x2 2x3 8.3 x1 x2 5x3 4.2 表 4.2.1 方程组进行运行结果 迭代方法 次数 精确解 近似解 Jacobi 迭代法 7 1.08030303030303 1.08004290000000 1.06493506493507 1.06474841000000 1.26904761904762 1.26874

47、770000000 Gauss-Seide 迭代法 5 1.08030303030303 1.08027154136992 1.06493506493507 1.06492268041955 1.26904761904762 1.26903884435789 SOR 方法 5 1.08030303030303 1.08027154136992 1.06493506493507 1.06492268041955 1.26904761904762 1.26903884435789 4.2.2 线性方程组运行结果分析 从运算的结果分析 SOR运行的结果与精确值相比较为准确，而且迭代的次数少。 在日常

48、生活中，迭代法常用的优 Jacobi 迭代法 ,Gauss-Seidel 迭代法， SOR方法 中 Jacobi 迭代法简单，并具有很好的串行算法，很适合并行计算，但收敛速度较慢。 Gauss-Seidel 迭代法是典型的串行算法，在 Jacobi 迭代法与 Gauss-Seidel 迭代法同 时收敛的条件下，后者比前者收敛的快，但两种迭代收敛收敛域互不相容，不能互相 替代。 SOR 方法是一种应用极为广泛的方法，但选取最佳松弛因子比较困难，常通过 试算来确定最佳松弛因子。各种方法都有其利弊，但是在计算方程组时要求精确度很 高时往往选用 SOR方法。 莹谐龌蕲賞组靄绉嚴减 。 16 附录 1、

49、复化梯形积分公式的源程序 function out1,out2,out3=fhtx(a,b,f) format long ; clc; if nargin=3 wc=0.5*10(-6); end % 控制输入参数结束 disp( fhtx? ? YD?) tic; n=1; s1=(subs(f, x ,a)+4*subs(f, 缚縟糶。 x ,(a+b)/2)+subs(f, x ,b)*(b-a)/n/6; 麸肃鹏镟轿騍镣 n=2; s2=s1+1; while abs(s2-s1)wc s1=s2; n=n*2; s=0; h=(b-a)/n; for i=1:n s=s+4*subs

50、(f, end ; for i=1:n-1 s=s+2*subs(f, end ; s=(subs(f, x s2=s; end ; time=toc; jsz=s; % 控制输出参数结束 if nargout=1 out1=jsz; end if nargout=2 out1=jsz; out2=n; end if nargout=3 out1=jsz; x ,a+(i-1/2)*(h); x ,a+i*h); ,a)+subs(f, x ,b)+s)*h/6; 17 out2=n; out3=time; end 2、复化辛普生源程序 function out1,out2,out3=xps(

51、a,b,f) format long ; clc; if nargin=3 wc=0.5*10(-6); end %此处填复化辛普生公式 tic; n=1;s1=(subs(f, x ,a)+4*subs(f, 鄖禎銣腻鰲锬。 x ,(a+b)/2)+subs(f, x ,b)*(b-a)/n/6; 納畴鳗吶 n=2;s2=s1+1; while abs(s2-s1)wc s1=s2; n=n*2; s=0; h=(b-a)/n; for i=1:n s=s+4*subs(f, end ; x ,a+(i-1/2)*(h); for i=1:n-1 s=s+2*subs(f, end ; s=

52、(subs(f, s2=s; end ; x ,a+i*h); x ,a)+subs(f, x ,b)+s)*h/6; jsz=s2; time=toc; % 控制输出参数开始 if nargout=1 out1=jsz; end if nargout=2 out1=jsz; out2=n; end if nargout=3 out1=jsz; out2=n; out3=time; end 18 3、复化牛顿科特斯源程序 function jsz,n,time=kts(a,b,f) format long ; clc; fd=zeros(1,5); wc=0.5*10(-3); %此处填复化科

53、特斯公式 tic; wqja=a;wqjb=b; for i=1:5 fd(i)=wqja+(i-1)*(wqjb-wqja)/4; end ; s1=(7*subs(f, x ,fd(1)+32*subs(f, fd(4)+7*subs(f, n=2; h1=(b-a)/n; x ,fd(5)*(wqjb-wqja)/90; ,fd(2)+12*subs(f, x ,fd(3)+32*subs(f, 風撵鲔貓铁频钙蓟纠庙。 s2=0; for i=1:n wqja=a+(i-1)*h1;wqjb=a+i*h1; for i=1:5 fd(i)=wqja+(i-1)*(wqjb-wqja)/4

54、; end s2=s2+(7*subs(f, x ,fd(4)+7*subs(f, end x ,fd(1)+32*subs(f, x ,fd(5)*(wqjb-wqja)/90; x ,fd(2)+12*subs(f, x ,fd(3)+32*subs(f, 灭嗳骇諗鋅猎輛觏馊藹。 while abs(s2-s1)wc s1=s2; n=n*2;h1=(b-a)/n; s2=0; for i=1:n wqja=a+(i-1)*h1;wqjb=a+i*h1; for i=1:5 fd(i)=wqja+(i-1)*(wqjb-wqja)/4; end s2=s2+(7*subs(f, x ,fd

55、(4)+7*subs(f, end x ,fd(1)+32*subs(f, x ,fd(5)*(wqjb-wqja)/90; x ,fd(2)+12*subs(f, x ,fd(3)+32*subs(f, 铹鸝饷飾镡閌赀诨癱骝。 end jsz=s2; time=toc; end 19 4、龙贝格积分公式的源程序 function t=Romberg(fname,a,b,e) %Romberg 法求函数的积分 %fname 是被积函数 ,a 是积分上限 ,b 是积分下限 ,e 为精度 tic if nargine i=i+1;h=h/2; T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h)*h/2; for j=1:i T(i+1,j+1)=4j*T(i+1,j)/(4j-1)-T(i,j)/(4j-1); end end %T t=T(i+1,j+1); toc 攙閿频嵘陣澇諗谴隴泸。 趕輾雏纨颗锊讨跃满賺。 5、 jacobi 迭代法求方程组