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

输出Ri+1 i=i+1 T2^(i+5)=T2^(i+4)/2+ 2j=1f(i+4bsinx

xdx，程序初始要求输出需要得到的精度，最后输出得到

输入精度eps，积分上下限b,a 定义函数f(x)= 1 x=0; sin(x)/x 其他。 计算Ti(0eps 2j−1 b−a 2i+5+1+𝑎) Si+4=4Ti+5/3-Ti+4/3 Ci+3=16Si+4/15-S3/15 Ri+2=16Ci+3/15-Ci+2/15

程序代码：

%定义函数

function f=sin1(x) if(x==0) f=1; else f=sin(x)/x; end

%龙贝格算法求数值积分 clc;clear all; flag=zeros(1,4); In=zeros(25,4);

eps1=input('what precision do you want to get?'); i=0;

%计算一些初始序列 while(flag(4)<2) for k=flag(1)+1:i+4 for j=1:2^(k-1)

In(k,1)=In(k,1)+sin1(a+(2*j-1)*(b-a)/(2^k)); end if k>1

In(k,1)=In(k,1)/2^k+In(k-1,1)/2; else

In(k,1)=(b-a)*(sin1(b)+sin1(a))/2; end end

flag(1)=i+4; for k=flag(2)+1:i+3

In(k,2)=4*In(k+1,1)/3-In(k,1)/3; end

flag(2)=i+3; for k=flag(3)+1:i+2

In(k,3)=16*In(k+1,2)/15-In(k,2)/15; end

flag(3)=i+2; for k=flag(4)+1:i+1

In(k,4)=*In(k+1,3)/63-In(k,3)/63; end

flag(4)=i+1; i=i+1; end

%根据递推式计算龙贝格序列 i=1;

while (abs(In(i+1,4)-In(i,4))>eps1)

for j=1:2^(i+4)

In(i+5,1)=In(i+5,1)+sin1(a+(2*j-1)/(2^(i+5))); end

In(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)=*In(i+3,3)/63-In(i+2,3)/63; i=i+1; end

%输出高精度计算值 vpa(In(i+1,4))

程序执行结果：