《积分方程的数值计算技巧》由会员分享，可在线阅读，更多相关《积分方程的数值计算技巧（8页珍藏版）》请在金锄头文库上搜索。
1、实验七 积分方程的数值计算方法1. 实验描述计算定积分的近似值，起始容差1.用组合梯形公式M=10计算。2.用组合辛普生公式M=5计算。3.用龙贝格积分计算。4.用自适应积分方法计算。2. 实验内容1. 用组合梯形公式M=10计算。Y N 子区间宽度输入区间开始结束图1. 组合梯形算法流程图将积分区间划分为宽度为的M个子区间，再将各区间的面积求和即可得到近似面积。2.用组合辛普生公式M=5计算。输出SY N 结束子区间宽度输入区间开始图2. 组合辛普森算法流程图将积分区间划分为宽度为的2M个子区间，再由辛普森公式将各区间的面积求和即可得到近似面积。3.用龙贝格积分计算。结束开始输出面积近似值输
2、入区间图3.龙贝格积分算法流程图.由递归梯形公式序列得到递归辛普森序列序列。.由递归辛普森序列序列得到递归布尔公式序列。.通过理查森改进提高误差项的阶数。4.用自适应积分方法计算结束将辛普森公式中的区间划分为两个相等的子区间开始图4.自适应积分算法流程图在辛普森公式基础上，将区间再进行划分，即为自适应积分。3. 实验结果及分析真实值S = 0.199714662161441. 用组合梯形公式M=10计算。面积近似值S1 =0.16965032127666。误差error1=0.03006434088479。2.用组合辛普森公式M=5计算面积近似值S2 =0.19966805529718。误差e
3、rror2=。3.用龙贝格积分计算。表1.龙贝格积分表J梯形公式辛普森公式布尔公式0-0.00199504693265 000001-0.02186444123413-0.02848757266796000020.016117544001270.028778205746400.0325959243073600030.152657188447420.198170403262800.209463216430560.212270633765850040.188002801390540.199784672371580.199892290312170.199740370849970.1996912325
4、6403050.196791377983860.199720903514970.199716652257860.199713864352240.199713760405190.19971378242654面积近似值S3 = 0.19971378242654。误差error3=。4.用自适应积分方法计算面积近似值S4 = 0.19971391278871。误差error4=。4. 结论自适应计分方法和龙贝格积分法都可以得到较高的精确度。组合梯形和组合辛普森可以增大M值，得到较高的精确度。附件（代码）1.f=inline(sin(4*x)*exp(-2*x),x);%定义函数A=0;t=0;S=0
5、;for k=1:9 t=k*0.3; A=A+0.3*feval(f,t);endS=A+0.15*(feval(f,0)+feval(f,3)%组合梯形的近似值RS=int(sin(4*x)*exp(-2*x),0,3)%真实值error=RS-S%真实值与近似值的误差2.f=inline(sin(4*x)*exp(-2*x),x);%定义函数S=0;X=zeros(1,11);h=0.3;for j=1:11 X(j)=(j-1)*h;endfor k=1:5 S=S+h/3*(feval(f,X(2*k-1)+4*feval(f,X(2*k)+feval(f,X(2*k+1);endR
6、S=int(sin(4*x)*exp(-2*x),0,3);%真实值Serror=RS-S%真实值与近似值的误差3.function R,quad,err,h=romber(f,a,b,n,tol)%INPUT-f is the intergrand input as astring f% -a and b are upper and lower limits of interation% -n is the maxium number of rows in the table% -tol is the tolerance%OUTPT-R is the Romberg table% -quad
7、 is the quadrature value% -err is the error estimate% -h is the smallest step size uesdM=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(feval(f,a)+feval(f,b)/2;while(errtol)&(Jn)|(J4) J=J+1; h=h/2; s=0; for p=1:M x=a+h*(2*p-1); s=s+feval(f,x); end R(J+1,1)=R(J,1)/2+h*s; M=2*M; for K=1:J R(J+1,K+1)=R(J+1,K
8、)+(R(J+1,K)-R(J,K)/(4K-1); end err=abs(R(J,J)-R(J+1,K+1);endquad=R(J+1,J+1);4.function Z=srule(f,a0,b0,tol0)%INPUT-f is the intergrand input as astring f% -a0 and b0 are upper and lower limits of intergration% -tol0 is the tolerance%OUTPUT-Z is a 1x6vector a0 b0 S S2 err tol1h=(b0-a0)/2;C=zeros(1,3)
9、;C=feval(f,a0 (a0+b0)/2 b0);S=h*(C(1)+4*C(2)+C(3)/3;S2=S;tol1=tol0;err=tol0;Z=a0 b0 S S2 err tol1; function SRmat,quad,err=adapt(f,a,b,tol)%INPUT-f is the intergrand input as astring f% -a and b are upper and lower limits of interation% -tol is the tolerance%OUTPT-SRmat is the Romberg table% -quad i
10、s the quadrature value% -err is the error estimate% -Initialize valuesSRmat=zeros(30,6);iterating=0;done=1;SRvec=zeros(1,6);SRvec=srule(f,a,b,tol);SRmat(1,1:6)=SRvec;m=1;state=iterating;while(state=iterating) n=m; for j=n:-1:1 p=j; SR0vec=SRmat(p,:); err=SR0vec(5); tol=SR0vec(6); if (tol=err) %Bisec
11、t interval,apply Simpsonrule %recursively,and determine error state=done; SR1vec=SR0vec; SR2vec=SR0vec; a=SR0vec(1); b=SR0vec(2); c=(a+b)/2; err=SR0vec(5); tol=SR0vec(6); tol2=tol/2; SR1vec=srule(f,a,c,tol2); SR2vec=srule(f,c,b,tol2); err=abs(SR0vec(3)-SR1vec(3)-SR2vec(3)/10; %Accuracy test if (errtol) SRmat(p,:)=SR0vec; SRmat(p,4)=SR1vec(3)+SR2vec(3); SRmat(p,5)=err; else SRmat(p+1:m+1,:)=SRmat(p:m,:); m=m+1; SRmat(p,:)=SR1vec; SRmat(p+1,:)=SR2vec; state=iterating; end end endendquad=sum(SRmat(:,4);err=sum(abs(SRmat(:,5);SRmat=SRmat(1:m,1:6);8 / 8
