龙贝格算法求数值积分程序说明及流程图

陆韶琦3110000441程序说明：本程序计算sinxxdxba，程序初始要求输出需要得到的精度，最后输出得到数值积分值。流程图：输入精度eps，积分上下限b,a计算Ti(0i5)T2^(i+1)=T2^i/2+f(2j−1b−a2i+1+𝑎)2^ij=1计算Si（0i4）Si=4Ti+1/3-Ti/3计算Ci(0i3)Ci=16Si+1/15-Si/15计算Ri（0i2）Ri=16Ci+1/15-Ci/15|Ri+1-Ri|epsT2^(i+5)=T2^(i+4)/2+f(2j−1b−a2i+5+1+𝑎)2i+4j=1Si+4=4Ti+5/3-Ti+4/3Ci+3=16Si+4/15-S3/15Ri+2=16Ci+3/15-Ci+2/15i=1i=i+1输出Ri+1false定义函数f(x)=1x=0;sin(x)/x其他。程序代码：%定义函数functionf=sin1(x)if(x==0)f=1;elsef=sin(x)/x;end%龙贝格算法求数值积分clc;clearall;flag=zeros(1,4);In=zeros(25,4);eps1=input('whatprecisiondoyouwanttoget?');i=0;%计算一些初始序列while(flag(4)2)fork=flag(1)+1:i+4forj=1:2^(k-1)In(k,1)=In(k,1)+sin1(a+(2*j-1)*(b-a)/(2^k));endifk1In(k,1)=In(k,1)/2^k+In(k-1,1)/2;elseIn(k,1)=(b-a)*(sin1(b)+sin1(a))/2;endendflag(1)=i+4;fork=flag(2)+1:i+3In(k,2)=4*In(k+1,1)/3-In(k,1)/3;endflag(2)=i+3;fork=flag(3)+1:i+2In(k,3)=16*In(k+1,2)/15-In(k,2)/15;endflag(3)=i+2;fork=flag(4)+1:i+1In(k,4)=64*In(k+1,3)/63-In(k,3)/63;endflag(4)=i+1;i=i+1;end%根据递推式计算龙贝格序列i=1;while(abs(In(i+1,4)-In(i,4))eps1)forj=1:2^(i+4)In(i+5,1)=In(i+5,1)+sin1(a+(2*j-1)/(2^(i+5)));endIn(i+5,1)=In(i+5,1)/2^(i+5)+In(i+4,1)/2;In(i+4,2)=4*In(i+5,1)/3-In(i+4,1)/3;In(i+3,3)=16*In(i+4,2)/15-In(i+3,2)/15;In(i+2,4)=64*In(i+3,3)/63-In(i+2,3)/63;i=i+1;end%输出高精度计算值vpa(In(i+1,4))程序执行结果：