龙贝格积分实验报告

二、Romberg积分法1.变步长Romberg积分法的原理复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要事先估计出部长。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。而用复化公式的截断误差来估计步长，其结果是步长往往过小，而且''()fx和(4)()fx在区间[,]ab上的上界M的估计是较为困难的。在实际计算中通常采用变步长的方法，即把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg积分法的思想。在步长的逐步分半过程中，要解决两个问题：1.在计算出NT后，如何计算2NT,即导出2NT和NT之间的递推公式；2.在计算出NT后，如何估计其误差，即算法的终止的准则是什么。首先推导梯形值的递推公式，在计算NT时，需要计算1N个点处的函数值在计算出NT后，在计算2NT时，需将每个子区间再做二等分，共新增N个节点。为了避免重复计算，计算2NT时，将已计算的1N个点的数值保留下来，只计算新增N个节点处的值。为此，把2NT表示成两部分之和，即2122211222111222111[()()2()]21[()()2()2((21))]21[()()2()]((21))]22NNNNkNNNNNkkNNNNNNkkThfafbfakhhfafbfakhfakhhfafbfakhhfakh由此得到梯形值递推公式22211((21))12NNNNNkTThfakh因此111212122,[()()],211,()22hhbaTfafbhhTThfah由复化梯形公式的截断误差有2''112''2222(),12(),12NNNNbaIThfabbaIThfab若''()fx变化不大时，即''''12()()ff，则有22241()2413NNNNNTTITTT式（2）表明，用2NT作为定积分I的近似值，其误差大致为21()3NNTT，因此其终止条件为2NNTT其中是预先给定的精度。2.Romberg积分公式将上述方法不断推广下去，可以得到一个求积分的序列，而且这个序列很快收敛到所求的定积分。记(0)NNTT，将区间N等分的梯形值。(1)NNTS，将区间N等分的Simpson(2)NNTC，将区间N等分的Cotes。(3)NNTR，将区间N等分的Romberg。由其可构造一个序列(){}kNT，次序列称为Romberg序列，并满足如下递推关系：(0)(0)(0)1211[()()],((21)),2222NNNkbababaTfafbTTfakNN(1)(1)24,1,2,41kkkkNNNkTTTk以上递推公式就是Romberg积分递推公式。3.Romberg积分程序1.置1N，精度要求，1hba；2.计算(0)1[()()]2baTfafb；3.置22NNhh，并计算(0)(0)211((21))222NNNkbabaTTfakNN；4.置,2,1;MNNNK5.计算(1)(1)2441kkkkMMMkTTT；6.若1M，则转（7）；否则置2MM，1kk转（5）；7.若()(1)11kkTT，则停止计算（输出()1kT），否则转（3）。4.Romberg积分法的应用function[T,n]=romb(f,a,b,eps)doubleR;ifnargin4,eps=1e-8;endh=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b));n=1;J=0;err=1;while(erreps)J=J+1;h=h/2;S=0;fori=1:nx=a+h*(2*i-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;fork=1:JR(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1);enderr=abs(R(J+1,J+1)-R(J+1,J));n=2*n;endR;T=R(J+1,J+1)End其中输入项：f为被积函数，ab为积分区间的端点值，ep为积分精度；输出项：T是逐次积分表值，n是迭代次数，R是最后积分值。4.1程序调用可以将被积分函数编成函数文件，也可以直接使用内联函数来表示被积分函数，示例如下:f=inline('1/(1+x.^2)','x');[T,n,R]=romb(f,2,9,1e-9)运行后得出其迭代次数，最终积分结果以及龙贝格积分矩阵如表2-1所示，迭代次数N=64，最终的积分值R=0.3530.表2-1龙贝格积分矩阵3.课本例题求解2111200001ln(1)ln(1)sin()(1),(2),(3),(4)1+1xxxdxdxdxdxxxxx1当迭代精度ep=1e-9的条件下，迭代次数N=32，迭代结果R=0.6931表2-2式1对应的龙贝格积分矩阵0.75000.00000.00000.00000.00000.00000.70830.69440.00000.00000.00000.00000.69700.69330.69320.00000.00000.00000.69410.69320.69310.69310.00000.00000.69340.69310.69310.69310.69310.00000.69320.69310.69310.69310.69310.69312当迭代精度ep=1e-9的条件下，迭代次数N=32，迭代结果R=0.2722.表2-3式2对应的龙贝格积分矩阵0.17330.00000.00000.00000.00000.00000.24880.27400.00000.00000.00000.00000.26650.27230.27220.00000.00000.00000.27080.27220.27220.27220.00000.00000.27180.27220.27220.27220.27220.00000.27210.27220.27220.27220.27220.27223对于积分10ln(1)xdxx，由于积分下限0为其奇点，理论上无法进行数值积分，本题中近似取下限为1*10-9来进行计算。当迭代精度ep=1e-9的条件下，迭代次数N=16，迭代结果R=0.2722.0.74270.00000.00000.00000.00000.00000.00000.48330.39690.00000.00000.00000.00000.00000.39050.35960.35710.00000.00000.00000.00000.36280.35360.35320.35320.00000.00000.00000.35550.35300.35300.35300.35300.00000.00000.35360.35300.35300.35300.35300.35300.00000.35310.35300.35300.35300.35300.35300.3530表2-4式3对应的龙贝格积分矩阵0.84660.00000.00000.00000.00000.82880.82280.00000.00000.00000.82410.82250.82250.00000.00000.82290.82250.82250.82250.00000.82260.82250.82250.82250.82254.对于积分20sin()xdxx，同样积分下限0为积分函数的奇点，理论上无法进行数值积分运算，本题中仍取积分下限近似为1*10-9进行计算。当迭代精度ep=1e-9的条件下，迭代次数N=16，迭代结果R=1.3708.表2-5式4对应的龙贝格积分矩阵1.28540.00000.00000.00000.00001.34981.37130.00000.00000.00001.36551.37081.37080.00000.00001.36951.37081.37081.37080.00001.37041.37081.37081.37081.3708