数值计算方法第五章第二节 复化求积公式

第二节复化求积公式复化求积公式的基本思想：

将区间[a,b]

分为若干个小子区间，在每个小子区间上使用低阶的Newton-Cotes公式。然后把它们加起来，作为整个区间上的求积公式。

一、复化求积公式

1、复化梯形公式复化梯形公式为截断误差分析：2、复化Simpson公式复化Simpson公式的截断误差为复化Simpson公式为Example1ApproximatetheintegralUsingtheCompositeTrapezoidalruleandCompositeSimpson’srulex=eps:0.01:pi;y=sin(x)./x;plot(x,y);legend('f(x)=sin(x)/x');functionrs=trapezoid(f,a,b,n)h=(b-a)/n;r=(feval(f,a)+feval(f,b))/2;forj=1:n-1x=a+j*h;r=r+feval(f,x);endrs=r*h;

数值试验复化梯形公式Matlab程序将此程序存于work目录中数值试验

n=8,tp=1.84784230644461n=16,tp=1.85091414036536n=32,tp=1.85168137241373n=64,tp=1.85187313510989n=128,tp=1.85192107295303n=256,tp=1.85193305723690n=512,tp=1.85193605329681n=1024,tp=1.85193680231110n=2048,tp=1.85193698956462n=4096,tp=1.85193703637800n=8192,tp=1.85193704808136复化梯形公式数值实验结果functionrs=simpson(s,a,b,n)h=(b-a)/n;r=feval(s,a)+feval(s,b);

forj=1:2:n-1x=a+j*h;r=r+4*feval(s,x);

endforj=2:2:n-2x=a+j*h;r=r+2*feval(s,x);

end

rs=r*h/3;复化Simpson公式Matlab程序

数值试验

将此程序存于work目录中>>formatlong>>f=inline('sin(x)/x');>>a=eps;b=pi;n=8;

%eps是Matlab最小正数

>>sp=simpson(f,a,b,n)在Matlab命令窗口中，进行如下操作：

数值试验数值试验n=8,sp=1.85195370083255n=16,sp=1.85193808500560n=32,sp=1.85193711642985n=64,sp=1.85193705600861n=128,sp=1.85193705223407n=256,sp=1.85193705199819n=512,sp=1.85193705198345n=1024,sp=1.85193705198253n=2048,sp=1.85193705198247n=4096,sp=1.85193705198246n=8192,sp=1.85193705198247复化Simpson公式数值实验结果

数值试验结果用复化Simpson公式求解，当区间分成64份时，计算结果：sp=1.85193705600861比用复化梯形公式，当区间分成8192份时的计算结果：tp=1.85193704808136还要精确。启发：寻求精度更高的求解方法！利用复合梯形公式、复合simpson公式、复合Cotes公式等计算定积分时，需根据精度要求确定节点数，即步长的大小，需要笔算求高阶导数解决办法：采用变步长算法通常采取将区间不断对分的方法，即取n=2k

，反复使用复合求积公式，直到相邻两次计算结果之差的绝对值小于指定的精度为止。二、变步长复化求积公式变步长复化求积公式的基本思想：将区间[a,b]逐次分半，建立递推公式，按递推公式计算，直到满足精度要求。1.变步长复化梯形公式下面推导由n到2n的复化梯形公式变步长复化梯形公式的递推公式：（由n到2n）变步长复化梯形求积公式的算法2.变步长复化Simpson公式同理可得变步长复化柯特斯公式二、龙贝格(Romberg)方法龙贝格(Romberg)算法是将理查逊(Richardson)外推法应用于数值积分，由低精度求积公式推出高精度求积公式的算法。龙贝格算法的计算公式k01234误差龙贝格序列计算流程龙贝格算法的计算过程逐次外推计算流程如下：k0123function

R=romberg(s,a,b,n)h=b-a;R(1,1)=h*(feval(s,a)+feval(s,b))/2;fori=2:nx=linspace(a,b,2^(i-1)+1);sm=0;fork=2:2:2^(i-1)sm=sm+feval(s,x(k));

endR(i,1)=R(i-1,1)/2+h*sm/2;forj=2:iR(i,j)=R(i,j-1)+(R(i,j-1)-R(i-1,j-1))/(4^(j-1)-1);

endh=h/2;end

%Romberg.mRomberg算法例2使用Romberg积分法计算>>formatlong>>f=inline('sin(x)/x');>>a=0.0+eps;b=pi;n=4;>>R=romberg(f,a,b,n)R=

1.570796326794900001.785398163397451.85693210893163001.835508123280791.852211443241901.85189673