科学计算与数学建模第三章.ppt

科学计算与数学建模,中南大学数学科学与计算技术学院,湘江流量估计模型,第三章湘江流量估计模型,3.1湘江水流量估计的实际意义,水流量是水文特征值的一个重要指标，而水文特征值对于水资源的合理利用，防洪以及抗旱具有指导性的作用，因此湘江水流量估计对于湘江流域的社会经济和人民生活具有重大的影响。现根据实际测量得到湘江某处河宽700m，其横截面不同位置某一时刻的水深如表3.1.1所示。若此刻湘江的流速为0.5m/s，试估计湘江此刻的流量。要计算湘江水流量就需要知道其横截面面积，如果知道此处江的水深曲线函数，则其横截面面积为。但是在实际中是不可能精确得到的，那么怎样求出足够高精度的横截面面积的近似值。,表3.1.1湘江某处横截面不同位置的水深数据单位：m,3.1.1数值求积的必要性在高等代数中，曾用牛顿-莱布尼兹（Newton-Leibniz）公式：（其中F(X)是f(x)的一个原函数）来计算定积分。但是，在工程技术和科学研究中，常常遇到如下情况：f(x)的结构复杂，求原函数困难；f(x)的原函数不能用初等函数表示；f(x)的精确表达式不知道，只给出了一张由实验提供的函数表。对于这些情况，要计算积分的精确值都是十分困难的，这就要求建立积分的近似计算方法。此外，积分的近似计算又为其它一些数值计算，例如微分方程数值解、积分方程数值解等提供了必须的基础。,其中为n次插值基函数。用近似代替被积函数f(x)，则得：,3.1.2构造数值求积公式的基本方法可以从不同的角度出发，通过各种途径来构造数值求积公式。但常用的一个方法是，利用插值多项式来构造数值求积公式。具体做法如下：在积分区间a,b上取一组点：作f(x)的n次插值多项式：,（3.1.1）,若记,得数值求积公式：,这样的求积公式称为机械求积公式。,（3.1.2）,（3.1.3）,其中称为求积节点，称为求积系数。若求积公式（3.1.3）中的求积系数是由（3.1.2）确定的，则称该求积公式为插值型求积公式。本章主要讨论插值型求积公式。,3.1.3求积公式的余项积分的真值与由某求积公式给出的近似之差，称为该求积公式的余项，记作。例3.1.1求积公式（3.1.3）的余项为：,如果求积公式（3.1.3）是插值型的，则由上知：,（3.1.4）,为了使一个求积公式能对更多的积分具有良好的实际计算意义，就应该要求它对尽可能多的被积函数都准确地成立。在计算方法中，常用代数精度这个概念来描述。定义3.1.1若求积公式：对任意不高于m次的代数多项式都准确成立，而对于却不能准确成立，则称该公式的代数精度为m。,例3.1.2梯形公式的代数精度m=1。一个求积公式的代数精度越高，它就越能对更多的被积函数f（x）准确（或较准确）地成立，从而具有更好的实际计算意义。由插值型求积公式的余项（3.1.4）易得：定理3.1.1含有n+1个节点的插值型求积公式（3.1.3）的代数精度至少为n。,（3.1.5）,3.1.4求积公式的代数精度,3.2数值求积的Newton-Cotes（牛顿-柯特斯）方法,在3.1中，介绍了插值型求积公式及其构造方法。在实际应用时，考虑到计算的方便，常将积分区间等分之，并取分点为求积分节点。这样构造出来的插值型求积公式就称为牛顿-柯特斯（Newton-Cotes）公式。本节在介绍一般牛顿-柯特斯公式的基础上，介绍几个常用的牛顿-柯特斯公式以及这些公式在实际计算时的用法。,3.2.1Newton-Cotes（牛顿-柯特斯）公式若将积分区间a,bn等分，取分点作为求积节点，并作变量替换,那么插值型求积公式(3.1.3)的系数由（3.1.2）可得：,记,（3.2.1）,则于是，由（3.1.3）就可写出相应的插值型求积公式：,（3.2.2）,这就是一般的牛顿柯特斯公式，其中称为柯特斯系数。从柯特斯系数的算式（3.2.1）可以看出，其值与积分区间a,b及被积函数f(x)都无关，只要给出了积分区间的等分数n，就能毫无困难地算出为了便于应用，部分柯特斯系数列见表3-1,例3.2.2当n=1时有两点公式：当n=2时有三点公式当n=4时有五点公式：,（3.2.4）,（3.2.5）,其中,求积公式（3.2.3）就是梯形公式。求积公式（3.2.4）称为辛普生（Simpson）公式其几何意义就是通过三点的抛物线围成的曲边梯形面积近似地代替原曲边梯形面积（见图3.2.1）。因此，求积公式（3.2.4）又名抛物线公式。求积公式（3.2.5）称为柯特斯公式。,梯形公式、辛普生公式和柯特斯公式，是三个最基本、最常用的等距节点下的求积公式。,下述定理给出了这些求积公式的余项。定理3.2.1若在a,b上连续，则梯形公式（3.2.3）的余项为：若在上连续，则辛普生公式（3.2.4）的余项为：若在上连续，则柯特斯公式(3.2.5)的余项为：其中,（3.2.6）,（3.2.7）,（3.2.8）,3.2.2复合Newton-Cotes（牛顿-柯特斯）公式,由定理3.2.1知，当积分区间较大时，直接使用牛顿-柯特斯公式所得积分近似值的精度是很难得到保证的。因此在实际应用中，为了既能提高结果的精度，又使算法简便且易在电子计算机上实现，往往采用复合求积的方法。所谓复合求积，就是先将积分区间分成几个小区间，并在每个小区间上用低阶牛顿柯特斯公式计算积分的近似值，然后对这些近似值求和，从而得到所求积分的近似值。由此得到的一些具有更大实用价值的数值求积公式，统称为复合求积公式。,例3.2.3先将区间a,bn等分，记分点为，其中，称为步长，然后在每个小区间上应用梯形公式（3.2.3），即：就可以导出复合梯形公式：,若将所得积分近似值记成，并注意到，则上式即为：仿上，可得复合辛普生公式：和复合柯特斯公式：,（3.2.9）,（3.2.10）,其中定理3.2.2若f(x)在积分区间a,b上分别具有二阶、四阶和六阶连续导数,则符合求积公式(3.2.9)、(3.2.10)和(3.2.11)的余项分别为：,（3.2.11）,其中，且当h充分小时，又有：,（3.2.14）,(3.2.12),（3.2.13）,(3.2.15),（3.2.16）,（3.2.17）,定义3.2.1对于复合求积公式，若当h0时有则称是p阶收敛的。定理3.2.3复合求积公式（3.2.9）、（3.2.10）和（3.2.11）分别具有二阶、四阶和六阶收敛性。证明由收敛性的定义，从（3.2.19）可以看出，复合梯形公式（3.2.9）具有二阶收敛性。同样，可证明复合辛普生公式（3.2.10）和复合柯斯特公式（3.2.11）分别具有四阶和六阶收敛性。,3.2.3误差的事后估计与步长的自动选择虽然可用余项公式（2.12）（2.17）来估计近似值的误差，也可以根据精度要求用这些公式来确定积分区间的等分数，即确定步长h。但由于余项公式中包含被积函数f(x)的高阶导数，在具体计算时往往会遇到困难。因此，在实际应用时，常常利用误差的事后估计法来估计近似值的误差或步长h。该方法的大致做法是：将积分区间逐次分半，每分一次就用同一复合求积公式算出相应的积分近似值，并用前后两次计算结果来判断误差的大小。其原理和具体做法是：对于复合梯形公式（3.2.9），由余项公式（3.2.12）或（3.2.15）可以看出，当f(x)在积分区间上变化不大或积分区间a,b的等分数n较大（即步长h较小）时，若将a,b的等分数改为2n（即将步长缩小到原步长h的一半），则新近似值T2n的余项约为原近似值余项的1/4，即：,其中I表示积分的真值。对I求解得：此式表明，若用T2n作为积分真值I的近似值，则其误差约为,(3.2.20),故在将区间逐次分半进行计算的过程中，可以用前后计算结果和来估计误差与确定步长。先算出Tn和T2n，若（为计算结果的允许误差），则停止计算，并取T2n作为积分的近似值；否则，将区间再次分半后算出新近似值T4n，并检查不等式是否成立，直到得到满足精度要求的结果为止。对于复合辛普生公式（3.2.10）和复合柯特斯公式（3.2.11），当所涉及的高阶导数在积分区间上,变化不大或积分区间的等分数较大时，由相应的余项公式可以看出：分别对I求解得,(3.2.21),(3.2.22),因此，估计新近似值和的误差，并判断计算过程是否需要继续进行，也可以像使用复合梯形法求积分近似值那样，在将积分区间逐次分半进行计算的过程中下去.3.2.4复合梯形法的递推算式上段介绍的变步长的计算方案，虽然提供了估计误差与选取步长的简便方法，但还没有考虑到避免在同一节点上重复计算函数值的问题，故有进一步改进的余地。先看复合梯形公式。,在利用（3.2.9）计算Tn时，需要计算n+1个点（它们是积分区间a,bn等分点的分点，不妨简称为“n分点”）上的函数值。当Tn不满足精度要求时，根据上面提供的计算方案，就应将各个小区间分半，计算出新近似值。若利用（3.2.9）进行计算T2n，就需要求出2n+1个点（它们是“2n分点”）上的函数值。而实际上，在这2n+1个2n分点中，包含有n+1个n分点，对应的函数值在计算Tn时早已算出。为了避免这种重复计算，下面分析近似值T2n与原有近似值Tn之间的联系。,由复合梯形公（3.2.9）知式：若注意到在2n分点：中，当k取偶数时是n分点，当k取奇数时，才是新增加的分点。将新增加的分点处的函数值从求和记号中分离出来，就有：,即,（3.2.23）,由递推公式（3.2.23）可以看出，在已经算出Tn的基础上再计算T2n时，只要计算n个新分点上的函数值就行了。与直接利用复合梯形公式（3.2.9）求T2n相比较，计算工作量几乎节省了一半。,3.3Romberg（龙贝格）算法,龙贝格（Romberg)算法是在积分区间逐次分半的过程中，对用复合梯形产生的近似值进行加权平均，以获得准确程度较高的一种方法，具有公式简练、使用方便、结果较可靠等优点。本节介绍它的基本原理和应用方法。,3.3.1Romberg（龙贝格）算法的基本原理上节中介绍的递推公式（3.2.23）或（3.2.24），虽然具有结构简单，易在电子计算机上实现等优点，但是由它产生的梯形序列，其收敛速度却是非常缓慢的。例3.3.1用此法计算的近似值时，要一直算到才获得误差不超过的近似值（见例3.2.5）。因此，用这种方法计算更复杂的高精度要求的积分近似值显然是费时、费力甚至是不可能的。,如何提高收敛速度，以节约计算工作量，自然是人们极为关心的课题。由近似等式（3.2.20），用T2n作为积分真值I的近似值，其误差约为。因此，如果用作为T2n的一种补偿，可以期望所得到的新近似值：有可能比更好地接近于积分的真值。,如在例3.2.5中，和是两个精度很差的近似值，但如果将它们按（3.3.1）作线性组合，所得到的近似值：却具有七位有效数字，其准确程度比还要高，而计算只涉及求九个点上的函数值，其计算工作量仅为计算的。那么，按（3.3.1）式作线性组合得到的新近似值，其实质又是什么呢？通过直接验证，易知，亦即：,（3.3.2）,这表明在收敛速度缓慢的梯形序列基础上，若将与按（3.3.2）作线性组合，就可产生收敛速度较快的辛普生序列：同理，从近似等式（3.2.21）出发，通过类似的分析，可以得到：故在辛普生序列的基础上，将与按（3.3.3）作线性组合，就可产生收敛速度更快的柯特斯序列：这种加速过程还可以继续下去。,（3.3.3）,例3.3.2通过与的线性组合，可以在柯特斯序列的基础上，产生一个称为龙贝格序列的新序列，即：经过进一步的分析，可以证明，当满足一定条件时，龙贝格序列比柯特斯序列更快地收敛到积分的真值。综上可知，可以在积分区间逐次分半的过程中利用公式（3.3.2）、（3.3.3）和（3.3.4），将粗糙的近似值逐步“加工”成越来越精确的近似值。也就是说，将收敛速度缓慢的梯形序列逐步地“加工”成收敛速度越来越快的新序列：。这种加速的方法就称为龙贝格算法。其加工过程如图3.3.1，其中圆圈中的号码表示计算顺序。,（3.3.4）,图3.3.1,3.3.2Romberg（龙贝格）算法计算公式的简化为了便于上机计算，引用记号来表示各近似值，其中k仍代表积分区间的二分次数，而下标m则指出了近似值所在序列的性质：当m=0时，在梯形序列中；当m=1时，在辛普生序列中；当m=2时在柯特斯序列中.,例3.3.4表3.3.1中的各近似值，若用记号表示，则如表3.3.2所示。,引入上面的记号后，龙贝格算法所用到的各个计算公式可以统一为：,（3.3.5）,相应的程序框图见图3.3.2。其中为最大二分次数。最后指出下列几点：（1）当较大时，由（3.3.5）第三式知。因此，在实际计算中，常规定，即在计算到出现龙贝格序列为止。在这种情况下，程序框图3.3.2应做相应的修改，需将“按式（3.3.5）计算”改为“按式（3.3.5）计算”，并将精度控法改为（2）为防止假收敛，可设置最小二分次数。当时，跳过精度判别而继续运算；（3）可以用二维数组来存放与参加运算，也可用一维数组。,图3.3.2,3.4Gauss(高斯)型求积公式与测量位置的优化选取,下面介绍一种高精度的求积公式高斯（Gauss）型求积公式。3.4.1Gauss（高斯）型求积公式的定义在3.2中，限定把积分区间的等分点作为求积节点，从而构造出一类特殊的插值型求积公式，即牛顿-柯特斯公式。这种做法虽然简化了计算，但却降低了所得公式的代数精度。,例3.4.1在构造形如的两点公式时，如果限定求积节点，那么所得插值型求值公式：的代数精度仅为1。,(3.4.1),(3.4.2),但是，如果我们对（3.4.1）中的系数A0,A1和节点X0,X1都不加限制，那么就可以适当选取A0,A1和X0,X1，使所得公式的代数精度m1。事实上，若只要求求积公式（3.4.1）对函数都准确成立，只要和满足方程组,解之得,（3.4.3）,代入（3.4.1）即得：容易验证，所得公式（3.4.4）是代数精度的插值型求积公式。同理，对于一般求积公式：,（3.4.4）,（3.4.5）,只要适当选择2n+2个待定参数和使它的代数精度达到2n+1也是完全可能的。定义3.4.1若形如（3.4.5）的求积公式代数精度达到了，则称它为高斯型求积公式，并称相应的求积节点为高斯点。,3.4.2Gauss（高斯）求积公式的构造与应用可以像构造两点高斯型求积公式（3.4.4）那样，通过解形如（3.4.3）的方程组来确定高斯点和，从而构造n+1点高斯型求积公式。但是，这种做法要解一个包含有2n+2个未知数的非线性方程组，其计算工量是想当大的。一个比较简单的方法是：,（1）先用区间a,b上的n+1次正交多项式确定高斯点（2）然后利用高斯点确定求积系数当积分区间是-1,1时，两点至五点高斯型求积公式的节点、系数T和余项见表3.4.1，其中利用表3.4.1，可以方便地写出相应的高斯型求积公式。例如，当n=2时，由表3.4.1知：,表3.4.1,故得两点高斯型求积公式：又如，当n=3时，由表3.4.1可以查出三个求积节点和对应的三个系数（注意，系数0.55555556应连用两次），从而得到三点高斯型求积公式：对于一般区间a,b上的积分，也可以利用表3.4.1写出高斯型求积公式。,先作变量替换，令将区间上的积分转化为区间-1,1上的积分：记则等式（3.4.6）右端的积分为,（3.4.6）,利用表3.4.1，对于给定的n=1,2,3,4，可以