Simpson-Rule-Summary---辛普森法则.doc

考虑积分，如果在区间a，b内取等间隔的N份，间隔长度为h，简述矩形（、梯形、Simpson法则）计算积分的i）理论、误差精度分析，和算法计算流程。解：对于缓变函数我们可以用各个区间中点上函数值作来近似该区间的平均值其中。矩形法则：f(x)在区间a,b上的积分用矩形求积定义如下第i个区间对积分的贡献为： 如果围绕该区间中点的邻域内对函数f(x)作泰勒级数展开， 有其中，和分别表示了f(x)在处的一阶，二阶和三阶导数。相应地，积分在子区间内的值可以表示为其中第一项是矩形积分的近似值，第二项则由于其中的积分等于零而消除。从而，矩形法则在宽度为h的单个子区间内的最高阶误差由第三项给出在整个a,b区间上的总误差则通过将所有N个子区间的贡献相加得到其中我们利用了Nh=(b-a)，并且取为f(x)在a,b上的二阶导数的均值。梯形法则： 如果用一个梯形来近似代替每个子区间的面积，梯形的四个顶点分别位于，与。梯形的面积为 （2-8）对于整个a,b区间，积分值由所有窄带的和给出 （2-9）与其他值相比，与的贡献只具有一半的权重，这是由于它们位于区间的两个端点。根据Euler-Maclaurin积分求和公式，我们有 （2-10）系数为伯努利数，其满足以下下规律：，且当时，。有 ，而其他所有奇数阶项为零。前几个偶数阶项的值如下 更一般的，它们可以通过来生成梯形法则的优点在f(x)的值只在格点（网点）上，即处求出。，此方法中由有限步长引起的误差的预期值是矩形法则中的2倍。而另一方面，梯形法则的简洁性则有助于进行进一步的改进。此方法可以与外插法一起使用而构成一种非常高效的数值积分算法。框图2-1列出了梯形法则算法的主要步骤。该算法还包含了一种方法，即由一个小的子区间数开始，通过循环迭代逐步增加数目，而无需在格点上重复计算被积函数的值。Box 2-1 Subroutine TRAPZ(A,B,KEY,RSLT,H)Trapezoidal rule integrationArguments: A: Lower integration limit. B: Upper integration limit. KEY=0: Initiate a new integration; 0: subsequent iterations. RSLT: Value of the integral. H: Integration step size.Initialization: Define the integrand as an external function.1. First call (KEY=0):(a) Calculate the integral using only two subintervals.(i) Define the number of subintervals to be 2.(ii) Let SUM be the average of the two end points plus the midpoint.(iii) Define the step size H to be half of the distance from A to B.(b) Change the value of KEY to 1.(c) Return the product of H and SUM.2. Subsequent calls (KEY0):(a) Calculate the contribution from points halfway between old ones:(i) Start from the middle of the first subinterval.(ii) Add the value of the integrand at this point to SUM.(iii) Move to the next subinterval and repeat steps (i) and (ii).(b) Double the number of subintervals and halve the step size H.(c) Return H times SUM.Simpson法则当积分的上下限相对于展开的中心点对称时，积分泰勒展开式中含有f(x)的奇数阶导数的项都将等于零。利用这一性质，我们可以在相邻的两个子区间内对面积作泰勒级数展开，可得积分的精确度可以达到。利用中心差分，可将f(x)的二阶导数近似表为将此结果代入前述中心查分表式，则子区间xi-1,xi+1内的积分近似为对于整个a,b区间，结果为假设区间的数目N为偶数，则奇数格点的贡献是偶数格点贡献的2倍。这种权重上的差异来自我们为修正基本方法所得的一级结果所引入的f(x)二阶导数的贡献。而两端点和的权重则仅为偶数点的一半。Simpson法则的另一优点在于其自然引出了一种算法，即通过迭代使积分达到所需要的精确度。误差估计-原则上，我们可以通过估计所省略的高阶项大小得出数值计算中的误差。其前提是知道高阶导数（如）的平均值），这样做并不容易。比较相继的两次迭代结果是较为简便的一种方法。例如，如果两次计算值Ia,b之差小于预设的允许值，则我们认为计算已经达到了所需的精确度。对于Simpson法则，该算法可通过如下方法实施。假设我们在x=a和x=b之间对函数f(x)进行积分，且要求达到精确度。各个节点上的函数值对积分的贡献可分为三个部分，即端点区，奇数点与偶数点的贡献，其中，。可用，及以如下形式写出积分由此可见，对于一个新的N，我们可利用上一次循环中已有的，及的值。为便于进行迭代计算，应从N=开始，其中应为某个较小而合理的区间数目，例如取=6。这是迭代计算的第一步。在下一轮迭代中，将区间的数目翻倍，而令步长h减半，即令N=2。以这种方式改变子区间数目的好处在于，对于一个新的N，来自两个端点的贡献（存储为）是不变的，而偶数点与上一次迭代计算过的偶数点完全相同。也就是说，在新的迭代中的的值就是上一次循环的与之和。所需计算的仅仅是函数在（新的）奇数点上的值。从而我们得到结果 其中。利用上面等式进行变步长的循环计算即。如果计算值与前一循环所得结果之差小于，则可认为计算结果已收敛；否则即再次令子区间数目加倍，使N=N3=2N2。然后依此类推，直到计算结果收敛或者达到最大迭代次数。框图给出了算法的主要步骤。Box 2-2 Program DM_SIMPSSimpsons rule integrationInitialization:(a) Set as error tolerance and N as the initial number of points.(b) Set a maximum for the number of iterations allowed.1. Input a as the lower limit of integration and b as the upper limit.2. Make an initial calculation of the integral with N points:(a) Calculate the contributions from end points.(b) Sum over the contributions from odd points.(c) Sum over the contributions from even points.(d) Calculate the integral from these three terms using Eq. (2-16). 3. Double the number of points and calculate the value of the integral again:(a) Sum the previous even and odd contributions as the new even value.(b) Calculate the contribution from the new odd points.(c) Calculate the integral with Eq. (2-16) using the new values.4. Compare the new result with the previous one.(a) If the difference is greater thango back to step 3.(b) If the difference is smaller or equal to, return the result.举一个简单的例子就可以说明Simpson法则的威力。考虑积分众所周知，其结果为cos(0)=1。为了用数值方法求出其值，我们由N=6开始迭代。除了分别位于0和/2的两个端点，我们还分别在=，及处有三个奇数点，在=，处有两个偶数点。在第一次循环中，求出了被积函数在这七个点上的值。表2-2给出了积分间隔h的大小及所得的积分值。Table 2-2: Example of the accuracy of Simpsons rule for different N.Nh61224481.00002621.00000150.99999990.99999990.2620.1310.0650.033exact value1.0000000在这个简单的例子中，我们实际上可以用（2-12）式来估计对任一给定的N的预期误差。在任一给定的子区间中，在Simpson法则已包含的二阶导数项之后的第一个非零项是对于整个积分，我们必须将此结果乘以N并求出f(x)在a,b区间内的四阶导数的平均值由于正弦函数的偶数阶导数仍然是正弦函数，我们得到|。在所讨论的例子中，(b-a)=