第2章-计算方法课件

第二章第二章 科学计算方法科学计算方法 科学研究、工程技术以及其它方面的大量的实际科学研究、工程技术以及其它方面的大量的实际应用都要涉及到应用都要涉及到计算计算。科学计算与科学实验，理论。科学计算与科学实验，理论研究共同成为科学方法论的基本环节研究共同成为科学方法论的基本环节.它们互相补充，它们互相补充，互相依赖，又相对独立互相依赖，又相对独立.其中，其中，数值方法的主要内容是研究求解数学模数值方法的主要内容是研究求解数学模型的算法及有关理论。型的算法及有关理论。科学计算过程描述为：科学计算过程描述为：实际实际问题问题数学数学模型模型数值数值方法方法计算结计算结果分析果分析 2.1 2.1 绪论绪论 计算机程序设计的核心是算法计算机程序设计的核心是算法.评价算法的两评价算法的两个主要标准是个主要标准是速度和精度速度和精度，速度涉及计算量，速度涉及计算量，精度涉及误差精度涉及误差.误差有误差有数据误差数据误差和和方法误差。方法误差。数据误差是由计数据误差是由计算机机器数字长有限所致；方法误差是由算机机器数字长有限所致；方法误差是由“近近似似”产生，近似是数值计算方法的一大特点。产生，近似是数值计算方法的一大特点。这是因数值计算时不必要有时也不可能取精确这是因数值计算时不必要有时也不可能取精确值。值。1.来源来源 一、一、误差的来源和分类误差的来源和分类 来源于来源于科学计算过程的四个主要的环节科学计算过程的四个主要的环节：数数学模型的建立、有关数据的获取、近似计算方学模型的建立、有关数据的获取、近似计算方法（数值方法）的确立、数值计算的进行。法（数值方法）的确立、数值计算的进行。2.分类分类 分为分为:模型误差、观测误差、截断误差和舍入误差模型误差、观测误差、截断误差和舍入误差 其中截断误差（又称方法误差）其中截断误差（又称方法误差）,它是由求其它是由求其近似近似 所致；舍入误差通常是四舍五入所致所致；舍入误差通常是四舍五入所致此两此两类误差即为数值分析所研究的计算误差类误差即为数值分析所研究的计算误差.二、绝对误差与相对误差二、绝对误差与相对误差 定义定义1 设实数设实数 x*为某一数据的准确值，它的为某一数据的准确值，它的近似值为近似值为 x，称称 e(x)=x x*为为 x 的的绝对误差绝对误差，简，简称误差称误差.绝对误差仅考虑近似值与准确值的差异，而相对误差绝对误差仅考虑近似值与准确值的差异，而相对误差还考虑数据本身的大小，相对误差常用百分数表示还考虑数据本身的大小，相对误差常用百分数表示.在实在实际测量或计算时，可根据具体情况估计际测量或计算时，可根据具体情况估计 e(x)或或 er(x)的大的大小范围小范围.当当 x*0 时，称时，称为为 x 的的相对误差相对误差.则称则称 (x)为为 x 的的绝对误差限绝对误差限.当当 x*0 时，称时，称定义定义2 设设 x 为准确值为准确值 x*的近似值，如果的近似值，如果 为为 x 相对误差限相对误差限.说明：说明：当绝对误差限当绝对误差限 (x)已知时，通常相对误差限已知时，通常相对误差限 r(x)仍然是未知的，所以在处理具体问题时，由于准仍然是未知的，所以在处理具体问题时，由于准确值未知，可以用近似值代替，即用确值未知，可以用近似值代替，即用 x 代替代替 x*。例例 1 零件的质量取决于参数零件的质量取决于参数 x，设设 x 的标定值为的标定值为 1 个个单位单位.生产中允许参数有一定误差，由此将零件分成生产中允许参数有一定误差，由此将零件分成 A，B，C 三个等级，等级由相对误差限决定，三个等级，等级由相对误差限决定，A等为等为1%，B 等为等为5%，C 等为等为10%.试确定三个等级的零试确定三个等级的零件参数允许变化的范围件参数允许变化的范围.A 等：等：x 0.99，1.01，B 等：等：x 0.95，1.05C 等：等：x 0.9，1.1.解解 由题设，标定值由题设，标定值 x*=1，根据相对误差限定义根据相对误差限定义将三个相对误差限分别代入，得将三个相对误差限分别代入，得有有 例例2 测量一物体的长度是测量一物体的长度是954cm，问测量数据的相对误差限问测量数据的相对误差限是多大？是多大？0.00053=0.053%解解 因为实际问题所截取的近似数，其因为实际问题所截取的近似数，其绝对误差限绝对误差限一般不超过最小刻度的半个单位一般不超过最小刻度的半个单位，所以当，所以当 x=954 cm 时，有时，有 (x)=0.5cm，而而x 的相对误差满足的相对误差满足 所以，测量数据的相对误差限可确定为所以，测量数据的相对误差限可确定为 r(x)=0.053%通通常常，若若近近似似值值 x 的的绝绝对对误误差差限限是是某某一一位位上上的的半半个个单单位位，该该位位到到 x 的的第第一一位位非非零零数数字字一一共共有有 n 位，则称近似值位，则称近似值 x 有有 n 位有效数字位有效数字.例如，考虑例如，考虑 的不同近似值的不同近似值.取三位数取三位数 x3=3.14，有有三、三、有效数字有效数字于是近似值于是近似值 x3 有有 3位有效数字位有效数字 例如例如 0.0031207，293.7048，两个数可表示为两个数可表示为 0.3120710 2，0.293704810 3 注注:浮点数在计算机程序中常用浮点数在计算机程序中常用 em 代替代替 10m。实数的实数的浮点表示法浮点表示法:0.a1a2an 10m 其中第一部分其中第一部分0.a1a2an是是尾数部尾数部 且a1非零非零,第第二部分二部分10m是是定位部定位部，用来确定小数点的位置用来确定小数点的位置.例如例如 0.31207e 2，0.2937048e+3表示规格化浮点数表示规格化浮点数 0.3120710 2，0.293704810 3 由定义由定义3易知：易知：若若 x 为准确值为准确值 x*的按的按“四舍五入四舍五入”原则而得到的近似值，且原则而得到的近似值，且 x 的第一位非零数字到末尾的第一位非零数字到末尾一共有一共有 n 位，则位，则 x 有有 n 位有效数字位有效数字.定义定义3 3 设设 x 可表示为可表示为规格化浮点数规格化浮点数形式形式 x=0.a1a2an 10m （2.1.1）其中，其中，a1，a2，an 都是都是0 9 中的任一整数，且中的任一整数，且 a1 0.若若 x 的绝对误差满足：的绝对误差满足：（2.1.2）则称近似数则称近似数 x 具有具有 n 位有效数字位有效数字.例如，例如，的近似值的近似值 x=3.142 有有4位有效数字位有效数字.定理定理 1 设设 x 是是有效数字位数为有效数字位数为 n 的近似数，其表的近似数，其表达式为式（达式为式（2.1.1）.则它的则它的相对误差满足相对误差满足 （2.1.3）反之，若近似数反之，若近似数 x 的相的相对误差满足对误差满足 （2.1.4）则则 x 至少至少有有n位有效数字位有效数字.x=0.a1a2an 10m （2.1.1）(2.13)定理定理1 揭示了数据的有效位数与该数据近似值的揭示了数据的有效位数与该数据近似值的相对误差之间关系相对误差之间关系.根据式（根据式（2.1.3），由于），由于a11，可知，一个具有可知，一个具有 n 位有效数字的近似数位有效数字的近似数 x，其相对，其相对误差满足误差满足 这说明，这说明，有效数位数与相对误差的阶码相对应有效数位数与相对误差的阶码相对应，x 的有效数字位数越多，相对误差越小的有效数字位数越多，相对误差越小.因此，要使相对误差限不超过因此，要使相对误差限不超过0.1%，只须不等式，只须不等式 10 n 0.001 成立即可成立即可.此时，取此时，取 n=3 即可即可.故当故当 x 取取 3 位有效数字位有效数字.例例3 要使要使 的近似值的近似值 x 的相对误差限不超过的相对误差限不超过0.1%，x 应取几位有效数字？应取几位有效数字？。解解:因因 的的第第一一位位数数字字为为5，所所以以近近似似数数 x 的的第第一一位位数字数字 a1=5，根据定理，根据定理1，当，当 x 取取 n 位有效数字时，有位有效数字时，有四、函数计算的误差估计四、函数计算的误差估计设设一一元元函函数数 y=f(x)，自自变变量量准准确确值值为为 x*，对对应应的的函函数数准准确确值值为为 y*=f(x*)，自自变变量量近近似似值值 x 的的误误差差为为 e(x)，误误差差限限为为 (x)，函函数数近近似似值值的的误误差差为为 e(y)，误误差限为差限为(y).利用一元函数全微分，得利用一元函数全微分，得由此得函数计算误差估计由此得函数计算误差估计当当 y0，x0 时，有相对误差估计时，有相对误差估计五、算术运算的误差估计五、算术运算的误差估计 设两个近似数设两个近似数 x1与与x2 的误差限分别为的误差限分别为 和和 ，则，则对这两个数的加、减、乘、除运算，可以利用对这两个数的加、减、乘、除运算，可以利用多元函数的误多元函数的误差估计差估计，得，得例例 4 设设 a=2.31，b=1.93，c=2.24 都是有三位有都是有三位有效数字的近似数，令效数字的近似数，令 p=a+bc，求，求 (p)和和 r(p)，并判断并判断 p 有几位有效数字有几位有效数字.解解 由于由于a，b，c都有三位有效数字，故都有三位有效数字，故 (a)=(b)=(c)=0.005 所以所以 (p)=(a)+(bc)(a)+b (c)+(b)c=0.005+1.930.005+2.240.005=0.02585 而而 p=2.31+1.932.24=6.6332 故故0.0039=0.39%因因(p)0.02585 0.05，故，故 p 中只有两位有效数字中只有两位有效数字.例例5 已测得某场地长和宽分别为已测得某场地长和宽分别为 x=110m，y=80m，已知已知 试求该场地面积的绝对试求该场地面积的绝对误差限和相对误差限误差限和相对误差限.相对误差限相对误差限解解:因面积因面积 S=xy，故绝对误差限，故绝对误差限=800.2+1100.1=27（m2）=0.0031=0.31%六、六、数值计算中的一些基本原则数值计算中的一些基本原则为使误差不增长，有一些基本原则为使误差不增长，有一些基本原则.1避免绝对值小的数作除数避免绝对值小的数作除数 这一原则主要指尽量避免除数绝对值远远小于被这一原则主要指尽量避免除数绝对值远远小于被除数绝对值的除法除数绝对值的除法.2避免两个相近的数据相减避免两个相近的数据相减 3要防止大数要防止大数“吃掉吃掉”小数小数 一个绝对值很大的数和一个绝对值很小的数直接相一个绝对值很大的数和一个绝对值很小的数直接相加时，很可能发生所谓加时，很可能发生所谓“大数吃小数大数吃小数”的现象，从而的现象，从而影响计算结果的可靠性影响计算结果的可靠性.这主要是计算机表示的数位这主要是计算机表示的数位数有限这一客观现实引起的数有限这一客观现实引起的.4尽量减少计算工作量尽量减少计算工作量 在考虑算法时应注意简化计算步骤，减少运算次数在考虑算法时应注意简化计算步骤，减少运算次数.计算机执行一个算法所花费的时间代价除了与问题的计算机执行一个算法所花费的时间代价除了与问题的规模大小有关外，规模大小有关外，主要依赖于计算过程中所用乘除法主要依赖于计算过程中所用乘除法次数的多少次数的多少，也就是计算工作量的大小，也就是计算工作量的大小.计算工作量小计算工作量小的算法不仅节约运行时间，而且使误差积累小的算法不仅节约运行时间，而且使误差积累小.5选用数值稳定性好的算法选用数值稳定性好的算法 对同一个数学问题，即使在数学公式已经确定的对同一个数学问题，即使在数学公式已经确定的情况下，仍然可以设计出不同的算法情况下，仍然可以设计出不同的算法.而不同的算法而不同的算法在执行过程中对数据误差的影响不一样，在执行过程中对数据误差的影响不一样，舍入误差舍入误差对计算结果影响不大的算法被称为对计算结果影响不大的算法被称为数值稳定的算法数值稳定的算法.方程：方程：f(x)=0如果如果 f(x)是非线性函数，则称为是非线性函数，则称为非线性方程非线性方程,其其中中f(x)是关于是关于x的一元函数的一元函数.2.2 非线性方程数值方法非线性方程数值方法 非线性方程的求解问题是科学与工程计算中常非线性方程的求解问题是科学与工程计算中常见的问题之一见的问题之一.但对高于但对高于4次的代数方程，不存在次的代数方程，不存在通用的求根公式，而对超越方程一般很难直接求通用的求根公式，而对超越方程一般很难直接求出其准确解出其准确解.所以，数值方法就是非常实用和有效的方法所以，数值方法就是非常实用和有效的方法.其其中，中，迭代迭代技术是一种常用技术技术是一种常用技术.其思想其思想是对根的是对根的逐逐次逼近次逼近.一般地，若用计算机求解非解线性方程有下列两步一般地，若用计算机求解非解线性方程有下列两步:第一步第一步 对方程的根进行隔离对方程的根进行隔离.找出找出隔根区间隔根区间（区间内只包含方程的一个根）（区间内只包含方程的一个根）.第二步第二步 利用利用迭代法计算迭代法计算满足一定精度的根近似值满足一定精度的根近似值.在方程的隔根区间内从给定的一个（或多个）出在方程的隔根区间内从给定的一个（或多个）出发值发值x0，按某种方法产生一个序列，按某种方法产生一个序列 x0，x1，x2，xn，.此序列在某种条件下收敛于方程的根此序列在某种条件下收敛于方程的根 x*.本节介绍非线性方程求根的本节介绍非线性方程求根的逐次逼近法逐次逼近法.同时也讨同时也讨论方法的论方法的收敛性收敛性和和误差估计误差估计等问题等问题.二分法本质上是一种二分法本质上是一种区间迭代算法区间迭代算法，在迭代过程中，在迭代过程中不断对隔根区间进行压缩，以区间中点逼近方程的根不断对隔根区间进行压缩，以区间中点逼近方程的根.它所涉及的理论是连续函数介值定理它所涉及的理论是连续函数介值定理:定理定理1 设函数设函数 f(x)在区间在区间a，b上连续，且上连续，且f(a)f(b)0，则，则f(x)=0 在区间在区间(a，b)内至少有一个根内至少有一个根.一一、二分法二分法方法方法(设函数设函数 f(x)满足定理满足定理1条件条件):1.求区间求区间a，b的中点的中点 x0=a+(b a)/2=(a+b)/22.计算计算f(x0),如果如果 f(x0)=0,则则 x0 是方程的解是方程的解;否则继续否则继续3.3.如果如果 f(x0)f(a)0，取，取 a1=a，b1=x0;如果如果 f(x0)f(b)0，取，取a1=x0，b1=b.此时区间此时区间 a1，b1的长度是的长度是a，b的一半的一半,将其作将其作为新的为新的a，b区间，重复区间，重复 1-3 直到找到根或找到满直到找到根或找到满足精度要求的近似根足精度要求的近似根(隔根区间充分小时隔根区间充分小时).说明说明 二分法计算过程中，设第二分法计算过程中，设第n区间为区间为an，bn，则，则(1)bn an=(b a)/2n (2)n充分大时，可取充分大时，可取xn=(an+bn)/2为方程的解为方程的解x*的近似值的近似值,且有如下且有如下收敛定理收敛定理 定理定理2 设设x*为方程为方程 f(x)=0在区间在区间a，b内的唯一根，内的唯一根，f(x)满足满足 f(a)f(b)0使对任意的使对任意的x0N(x*)，xn=(xn-1)产生的序列产生的序列 xn均均收敛于收敛于 x*，则称，则称迭代法局部收敛迭代法局部收敛.算法算法1（二分法二分法求解非线性方程算法）求解非线性方程算法）第四步：第四步：若若|b a|1 则转第二步；否则，输出则转第二步；否则，输出 x0 结束结束.第一步：第一步：输入误差限输入误差限 0，1，计算，计算 y1 f(a)，y2 f(b)；第二步：第二步：计算计算 x0 0.5(a+b)，y0f(x0)，若，若|y0|0，则，则输输 出出 x0，结束，结束.否则转第三步；否则转第三步；第三步：第三步：若若 y0 y1 0，则置，则置 b x0，y2 y0；否则；否则 a x0，y1 y0，转第四步；，转第四步；例例1 用二分法求下列方程在区间用二分法求下列方程在区间0，1内的一个根，要求内的一个根，要求 误差不超过误差不超过 1/25 所以所以 f(x)在区间在区间 0，1 上恰有一个根上恰有一个根.计算结果如下表计算结果如下表.解解:令令,计算得计算得n 函数值符号函数值符号 有根区间有根区间1 f(0.5)=0.1006 0 0.25，0.53 f(0.3755)=0.1317 0 0.375，0.54 f(0.4375)=0.0113 0 0.4375，0.5 取取 x4=0.5(0.4375+0.5)=0.4688 作为作为 x*的近似的近似值值.误差不超过误差不超过1/25.二、二、牛顿（牛顿（Newton）迭代法）迭代法算法思想算法思想 将方程将方程 f(x)=0 中函数中函数 f(x)在根的附近在根的附近线性化线性化，以线性方程的解逼近非线性方程的解以线性方程的解逼近非线性方程的解.理论依据理论依据 设函数设函数f(x)在有根区间在有根区间a，b二次连续可微，二次连续可微，x0是根是根 x*的一个近似值，则的一个近似值，则 f(x)在在 x0 处的泰勒展式为处的泰勒展式为 其中其中 在在 x 和和 x0之间之间.即即于是于是 f(x)=0 可近似转换可近似转换为为迭代公式的推导迭代公式的推导设设 f(x0)0，其解记为，其解记为这是较这是较 x0 更靠近更靠近x*的新的的新的近似值近似值.一般，在一般，在 xn 附近的附近的线性化方程为线性化方程为:设设 f(xn)0,其解记为其解记为n=0，1，Newton迭代公式迭代公式 这就是这就是,被称为牛顿迭代法的计算格式被称为牛顿迭代法的计算格式.选一个选一个初始点初始点 x0，由迭代公式便可得到迭代序列，由迭代公式便可得到迭代序列 x0，x1，xn，.例例3 用牛顿迭代法计算用牛顿迭代法计算 7 的平方根，记录并分析计的平方根，记录并分析计算过程中数据变化规律算过程中数据变化规律.解解 设设 ，得方程，得方程 x2 7=0.利用牛顿迭代利用牛顿迭代法，得计算格式法，得计算格式n=0，1，取初始值为取初始值为 x0=2.5，迭代计算的数据结果见下表，迭代计算的数据结果见下表 k xk xk xk 1 1 2.65000000000000 1.5000e-001 2 2.64575471698113 -4.2453e-003 3 2.64575131106678 -3.4059e-006 4 2.64575131106459 -2.1902e-012 5 2.64575131106459 0 由此可见，第五次迭代和第四次迭代数据的由此可见，第五次迭代和第四次迭代数据的15位数完位数完全相同，说明数据全相同，说明数据2.64575131106459 准确到小数点后准确到小数点后14位位.而得到这样的精度的数据只有而得到这样的精度的数据只有4次迭代次迭代.对于牛顿迭代法对于牛顿迭代法,迭代函数为迭代函数为 假定假定x*是是 f(x)=0 的单根，即的单根，即f(x*)=0，f(x*)0.由由 ，利用泰勒展式知，牛顿迭代法在根，利用泰勒展式知，牛顿迭代法在根x*附近至少平方附近至少平方收敛收敛.定理定理3 假设在的某邻域内具有连续的二阶导数，且设假设在的某邻域内具有连续的二阶导数，且设f(x*)=0，f(x*)0，则对充分靠近的初始值，则对充分靠近的初始值x0，牛顿，牛顿迭代法产生的序列迭代法产生的序列xn至少至少平方收敛平方收敛于于x*.牛顿迭代法的缺陷牛顿迭代法的缺陷:1.可能发生被零除错误可能发生被零除错误 当函数在它的零点附近，导函数的绝对值非常当函数在它的零点附近，导函数的绝对值非常 小时，运算会出现被零除错误；小时，运算会出现被零除错误；2.可能出现死循环可能出现死循环 当函数在它的零点有拐点时，可能会使迭代当函数在它的零点有拐点时，可能会使迭代 陷入死陷入死 循环循环.下面是两个牛顿迭代法不成功的下面是两个牛顿迭代法不成功的反例反例例例4 用牛顿迭代法解方程用牛顿迭代法解方程 x e x=0.同时可计算出当自变量增大时函数同时可计算出当自变量增大时函数 f(x)迅速趋近于零，迅速趋近于零，例如例如 f(x15)=0.0000000536.算法有可能错误地将算法有可能错误地将 x15做为做为根根.解解:令令 f(x)=xe x，f(x)恒为正，在区间恒为正，在区间1，内单内单调递减调递减.取取x0=2，用牛顿迭代法计算得，用牛顿迭代法计算得x1=4，x2=5.33333，x15=19.72354943，.发生被零除错误发生被零除错误例例5 用牛顿迭代法解方程用牛顿迭代法解方程 x3 x 3=0.解：解：对于对于 f(x)=x3 x 3，取，取x0=0，用牛顿迭代法计，用牛顿迭代法计 算，得算，得x1=3，x2=1.9615，x3=1.1472，x4=0.0066，.数列中的数列中的项以四项项以四项为一个周为一个周期重复期重复(见见右图右图)，算，算法将陷入法将陷入一种死循一种死循环中环中.定定理理4 4（非非局局部部收收敛敛性性定定理理）给给定定方方程程 f(x)=0，如如果函数果函数 f(x)C2a，b，且在，且在a，b满足条件：满足条件：（1）f(a)f(b)0.则则方方程程 f(x)=0 在在a，b 上上有有唯唯一一根根 x*，且且由由初初值值x0按按牛顿迭代公式牛顿迭代公式 n=0,1,2,求得的序列求得的序列 xn 二阶收敛二阶收敛于于x*.收敛性要求初始点要取得合适，否则会导致错误结果收敛性要求初始点要取得合适，否则会导致错误结果.定理定理4的解释的解释:1.定理定理4 中条件（中条件（1）保证了根的存在；）保证了根的存在；2.条件（条件（2）要求）要求 f(x)是单调函数，且是单调函数，且f(x)在在a，b上是上凸或下凸；上是上凸或下凸；3.条件（条件（3）取）取 x0a，b使使 f(x0)f”(x0)0 保证了保证了当当 xna，b时，有时，有 xn+1a，b.注记如下：注记如下：如如果果 f(x)的的二二阶阶导导数数大大于于零零，则则函函数数图图形形是是下下凸凸曲曲线线，根根据据条条件件（3），应应取取牛牛顿顿迭迭代代 的的 初初 始始 点点 使使 得得 f(x0)0（见右图）（见右图）y x O x0 x1 x2 如如果果 f(x)的的二二阶阶导导数数小小于于零零，则则函函数数图图形形是是上上凸凸曲曲线线，根根据据条条件件（3），应应取取牛牛顿顿迭迭代代的的初初始始点点使得使得 f(x0)0,当当 时时,由由弦弦截截法法产产生生的的序序列列 xn 收收敛敛于于 x*,且且收收敛敛阶至少为阶至少为1.618.2.3 线性方程组的直接解法线性方程组的直接解法 用用 E1,E2,En表示方程的行号表示方程的行号,则则 n 阶线性方程组可阶线性方程组可写成如下形式写成如下形式:写为矩阵形式写为矩阵形式 Ax=b (2.3.2)(2.3.1)分别称为方程组的分别称为方程组的系数矩阵系数矩阵、右端向量右端向量和和解向量解向量.线线性性方方程程组组的的数数值值解解法法分分为为直直直直接接接接法法法法和和迭迭代代法法两两类类.直直接接法法在在没没有有舍舍入入误误差差的的情情况况下下,通通过过有有限限步步计计算算求求得得方方程程组组的的准准确确解解.对对于于中中、小小型型方方程程组组(n1000)及及某某些些大大型型的的稀稀疏疏方方程程组组非非常常实实用用.实实际际计计算算中中,由由于于误误差差的的影影响响,直直接接法法只只能能求求得得方方程程组组的近似解的近似解.而而迭迭代代法法则则是是建建立立迭迭代代格格式式,从从初初始始向向量量出出发发,产产生生向向量量序序列列,逐逐次次逼逼近近方方程程组组的的准准确确解解.对对大大型型方程组方程组,常选用迭代法常选用迭代法.求解方程组求解方程组(2.3.1)的的 Gram 法则法则是公式求解方法是公式求解方法.对对一个一个n 阶方程组阶方程组,Gram 法则需要计算出包括法则需要计算出包括A的行列式的行列式在内的在内的(n+1)个个n 阶行列式阶行列式,而每一个行列式的计算需用而每一个行列式的计算需用n!(n-1)次乘法次乘法.对于对于n=20的中等规模问题的求解的中等规模问题的求解,也要耗也要耗费巨大机时费巨大机时;而作为一类最基本的直接法而作为一类最基本的直接法高斯消元高斯消元法法,其计算工作量只有大约其计算工作量只有大约2n3/3次浮点数运算次浮点数运算.一、一、高斯消元法高斯消元法基本思想基本思想 将一般的线性方程组将一般的线性方程组(2.3.1)经初等变换经初等变换化为等价的化为等价的上三角形方程组上三角形方程组,然后对上三角形方程然后对上三角形方程组直接求解组直接求解.1.上三角方程组求解的回代过程上三角方程组求解的回代过程如下形式的线性方程组称为上三角形方程组如下形式的线性方程组称为上三角形方程组 (2.3.3)当当 a11a22ann0 时时,求解过程为求解过程为 从最后一个方程开从最后一个方程开始始,逐个计算出未知元逐个计算出未知元(称为回代过程称为回代过程),即即算法算法1第一步：输入方程组系数矩阵和右端向量第一步：输入方程组系数矩阵和右端向量;回代过程所用的乘除法次数总共为回代过程所用的乘除法次数总共为 n(n+1)/2第二步：第二步：;1第三步：第三步：第四步：判断第四步：判断,若若k 1,则则k k-1,转第三步转第三步,否则否则,转第转第 五步五步;第五步：输出第五步：输出 ,结束结束.二、高斯消元过程二、高斯消元过程基本思想基本思想 对方程组的增广矩阵做行初等变换使其为对方程组的增广矩阵做行初等变换使其为上上三角形方程组三角形方程组,再再回代回代求出方程组的解求出方程组的解.例例1 用高斯消元法解方程组用高斯消元法解方程组:解解 Ei 表示增广矩表示增广矩阵的第的第 i 行行,用用 Ei+Ej Ej,表示第表示第 i 行乘数行乘数 加至第加至第 j 行行.该方程组的增广矩阵为该方程组的增广矩阵为:由回代求出方程组的解由回代求出方程组的解:x1=13 x2=8 x3=2已得到三角形方程组已得到三角形方程组 说明说明:对于方程组对于方程组(2.3.1),将其转化为等价的上三角将其转化为等价的上三角形方程组的过程称为形方程组的过程称为消元过程消元过程.算法算法2(消元过程的算法消元过程的算法)设方程组设方程组(2.3.1)的增广矩阵为的增广矩阵为第一步：输入第一步：输入 A=(aij)nn 和和 b=b1,bnT,k0;第二步：第二步：kk+1,i k;第三步：第三步：ii+1,计算计算 mik=aik/akk,bi bi mikbk,aij aij mik akj (j=k+1,n);如果如果 i n,则转第三步则转第三步,否则转第四步否则转第四步;第四步：如果第四步：如果k 1,转第五步第五步;否否则输出出 f1,f2,fn ,结束束.第四步：判断第四步：判断,如果如果 kn,转第三步转第三步;否则否则,jn,转第五步转第五步;随随着着计计算算技技术术的的发发展展,计计算算机机的的存存储储量量日日益益增增大大,计计算算速速度度也也迅迅速速提提高高,在在计计算算机机上上可可以以求求解解的的线线性性方方程程的的规规模模也也越越来来越越大大,在在实实际际应应用用中中,常常常常遇遇到到大大型型稀稀疏疏线线性性方方程程组组的的求求解解问问题题,因因此此寻寻求求能能保保持持稀稀疏疏性性的的有有效效算算法法就就成成为为数数值值代代数数中中的的一一个个非非常常重重要要的的课课题题.而而迭迭代代法法就就是其中之一是其中之一.2.4 线性方程组的迭代解法线性方程组的迭代解法迭代法思想迭代法思想:构造收敛的向量序列构造收敛的向量序列,使该序列收敛于所使该序列收敛于所求解方程组的准确解求解方程组的准确解.具体地具体地设设方方程程组组 Ax=b有有唯唯一一解解 x*,将将 Ax=b 变变形形为为等等价价的的方方程组程组 x=Bx+f 由此建立迭代公式由此建立迭代公式 x(k+1)=Bx(k)+f (k=0,1,2,).给定初始向量给定初始向量 x(0),按此公式计算的近似解向量序列按此公式计算的近似解向量序列 x(k),称此方法为称此方法为迭代法迭代法.若对任意初值若对任意初值 x(0),当迭代次当迭代次数无限增加时数无限增加时,序列序列 x(k)都有相同的极限都有相同的极限 x*,即即 x*=Bx*+f,则称迭代法是收敛的则称迭代法是收敛的,否则称为发散的否则称为发散的.迭代格式中的矩迭代格式中的矩阵阵B称为迭代矩阵称为迭代矩阵.一、雅可比迭代和高斯一、雅可比迭代和高斯-赛德尔迭代赛德尔迭代 迭代法适用于大型稀疏矩阵的线性方程组迭代法适用于大型稀疏矩阵的线性方程组,为了介绍该为了介绍该迭代法的思想迭代法的思想,考查下面三阶方程的迭代计算过程考查下面三阶方程的迭代计算过程.例例1 已知已知试用迭代法产生向量序列逐步逼近准确解试用迭代法产生向量序列逐步逼近准确解(注注:方程组的方程组的准确解为准确解为 x*=1,1,1T).解解 将原方程做等价变形将原方程做等价变形,由第一个方程解出由第一个方程解出 x1,由第二由第二个方程解出个方程解出 x2,由第三个方程解出由第三个方程解出 x3,得如下方程组得如下方程组:由此建立迭代格式由此建立迭代格式取初始向量取初始向量 x(0)=0,0,0T,迭代情况如表迭代情况如表2-6所示所示.说明说明:1.表表2-6中第一列数据为第一次迭代的结果中第一列数据为第一次迭代的结果,记记 为为 x(1);以以 x(1)为输入数据进行第二次迭代为输入数据进行第二次迭代,其结果为表中第其结果为表中第二列数据二列数据,记为记为 x(2);.当当 k=4时时,得近似解得近似解 x(4)=0.9987,0.9988,0.9991T.2.由上可知由上可知,对给定的方程组对给定的方程组,迭代法是利用迭代格式迭迭代法是利用迭代格式迭代计算出向量序列代计算出向量序列 x(1),x(2),x(n),逐步逼近方程组的准确解逐步逼近方程组的准确解.例例1所用迭代法称为所用迭代法称为雅可比雅可比(Jacobi)迭代法迭代法.k1234x1(k)0.77780.96300.99290.9987x2(k)0.80000.96440.99350.9988x3(k)0.86670.97190.99520.9991表表 2-6一般地一般地,n 阶线性代数方程阶线性代数方程 Ax=b(A R nn,b R n,x R n)可以写成如下形式可以写成如下形式:设设 aii 0 (i=1,2,n),由第由第i个方程解出个方程解出xi(i=1,2,n),得到与原方程组等价的方程组得到与原方程组等价的方程组(i=1,2,n).(i=1,2,n).将将迭代格式迭代格式(2.4.1)推广推广,有有(i=1,2,n).(2.4.2)从而得到雅可比迭代法的矩阵形式从而得到雅可比迭代法的矩阵形式x(k+1)=BJ x(k)+fJ ,k=1,2,n其中其中这就是雅可比迭代法的分量形式这就是雅可比迭代法的分量形式,即即将原方程组系数矩阵将原方程组系数矩阵A中主对角元分裂中主对角元分裂,设设 A=D L U 其中其中,D=diag(a11,a22,ann),aii0 (i=1,2,n),由于由于D 1 存在存在,于是于是BJ=D 1(L+U)=D 1(D A)=I D 1 A,fJ=D 1 b.称称 BJ =I D 1 A为为雅可比迭代矩阵雅可比迭代矩阵.算法算法1 (Jacobi 迭代算法迭代算法)(1)输入输入A,b,x(0),维数维数n,精度精度,最大容许迭代次数最大容许迭代次数N;(2)置置k=1;(3)计算计算 (i=1,2,n);(4)若若|x x(0)|,输出输出x,停机停机;否则转否则转(5);(5)若若 kN,则则 k k+1,xi(0)xi(i=1,2,n),转转(3);否则否则,输出失败信息输出失败信息,停机停机.高斯高斯-赛德尔赛德尔(Gauss-Seidel)迭代法迭代法引引入入 将将例例1中中的的方方程程组组及及迭迭代代格格式式(2.4.1)做做如如下下修修改改：尽尽量量利利用用最最新新计计算算数数据据,如如计计算算 时时,由由于于已已算算出出 ,所以用所以用 ,而不用而不用 ,则迭代格式则迭代格式(2.4.1)变为变为(2.4.3)取初始取初始值x(0)=0,0,0T,迭代情况如下表所示迭代情况如下表所示 k1234x1(k)0.77780.98390.99941.0000 x2(k)0.87780.99610.99981.0000 x3(k)0.97700.99870.99991.0000k=4时时,满足满足 =5.557810 4.说明说明:1.我们称迭代格式我们称迭代格式(2.4.3)为解为解Ax=b的的高斯高斯-赛德赛德尔尔(Gauss-Seidel)迭代法迭代法.2.对一般对一般n阶方程组阶方程组Ax=b,aii 0 (i=1,2,n),迭代格式迭代格式(2.4.3)就成为就成为这就是这就是Gauss-Seidel迭代法的分量形式迭代法的分量形式,其其矩阵形式矩阵形式为为x(k+1)=D 1(L x(k+1)+U x(k)+b)=D 1 L x(k+1)+D 1 U x(k)+D 1b.即即(i=1,2,n).(2.4.4)可改写为可改写为x(k+1)=(D L)1 U x(k)+(D L)1b =BG-S x(k)+fG-S (k=1,2,n)其中其中 BG-S=(D L)1 U 称为称为Guss-Seidel迭代矩阵迭代矩阵.计算实践表明计算实践表明,对许多问题对许多问题 Gauss-Seidel 迭代法迭代法确实比确实比 Jacobi 迭代法收敛快迭代法收敛快,但也有但也有Gauss-Seidel迭迭代比代比Jacobi迭代收敛慢的情况迭代收敛慢的情况,甚至还有甚至还有Jacobi迭代迭代收敛收敛,而而Gauss-Seidel迭代发散的情形迭代发散的情形.(1)输入输入A,b,x(0),维数维数n,精度精度,最大容许迭代次数最大容许迭代次数N;(2)置置k=1;算法算法2(Gauss-Seidel迭代算法迭代算法)(3)计算计算(i=2,n 1)(4)若若|x x(0)|,输出输出x,停机停机;否则转否则转(5);(5)若若kN,则置则置k k+1,xi xi(0)(i=1,2,n),转转(3);否则否则,输出失败信息输出失败信息,停机停机.定理定理1 1 迭代法迭代法(2.4.4)收敛的充分必要条件是迭代矩阵的收敛的充分必要条件是迭代矩阵的谱半径谱半径(B)1.收敛定理收敛定理 定理定理2 若若|B|1,则对于迭代格式则对于迭代格式 x(k+1)=B x(k)+f ,k=0,1,2,有有(2.4.5)(2.4.6)其中其中x*为方程组为方程组 Ax=b 的精确解的精确解.|为矩阵的算子范数为矩阵的算子范数,有有 定定理理3 若若A x=b的的系系数数矩矩阵阵A为为严严格格对对角角占占优优矩矩阵阵,则则Jacobi迭代和迭代和Gauss-Seidel迭代均收敛迭代均收敛.则称则称A为对角占优矩阵为对角占优矩阵.若上式中均为严格不等式若上式中均为严格不等式,则称则称A为严格对角占优矩阵为严格对角占优矩阵.定义定义2 设设 ,若若定理定理4 若若A为对称正定矩阵为对称正定矩阵,则则Gauss-Seidel迭代收敛迭代收敛.常微分方程数值解n1 引言引言n2 简单的简单的单步法单步法及基本概念及基本概念n3 Runge-Kutta方法方法1 引言2简单的单步法及基本概念简单的单步法及基本概念n2.1 Euler法，后退法，后退Euler法与梯形法法与梯形法n2.2 单步法的局部截断误差单步法的局部截断误差n2.3 改进改进Euler法法Euler法隐式隐式Euler法法若对初值问题积分形式采用右矩形公式可得：若对初值问题积分形式采用右矩形公式可得：称为隐式（后退）称为隐式（后退）Euler法。法。梯形法梯形法若对初值问题中对应的积分形式采用梯形公式可得：若对初值问题中对应的积分形式采用梯形公式可得：改进改进Euler法法将将梯形法和梯形法和Euler法相结合法相结合，可得到改进的，可得到改进的Euler法：法：3 Runge-Kutta方法n3.1 显式显式Runge-Kutta法法的一般形式的一般形式n3.2 二、三阶显式二、三阶显式R-K方法方法n3.3 四阶四阶R-K方法方法二阶二阶Runge-Kutta法：法：每推进一步须计算每推进一步须计算两两次函数值的次函数值的2阶单步法。阶单步法。三阶三阶Runge-Kutta法：法：每推进一步须计算每推进一步须计算三三次函数值的次函数值的3阶单步法。阶单步法。四阶四阶Runge-Kutta法：法：每推进一步须计算每推进一步须计算四次函数值的四次函数值的4阶单步法。阶单步法。