积分方程的数值计算技巧

实验七积分方程的数值计算方法
1. 实验描述
计算
32
sin(4)x
x e dx
-
⎰定积分的近似值，起始容差00.00001
ε=
1.用组合梯形公式M=10计算。

2.用组合辛普生公式M=5计算。

3.用龙贝格积分计算。

4.用自适应积分方法计算。

2. 实验内容
1. 用组合梯形公式M=10计算。

图1. 组合梯形算法流程图
将积分区间[],a b划分为宽度为()/
=-的M个子区间，再将各区间的面积
h b a M
求和即可得到近似面积。

2.用组合辛普生公式M=5计算。

Array
图2. 组合辛普森算法流程图
将积分区间[],a b划分为宽度为()/2
h b a M
=-的2M个子区间，再由辛普森公式将各区间的面积求和即可得到近似面积。

3.用龙贝格积分计算。

图3.龙贝格积分算法流程图
①.由递归梯形公式序列得到递归辛普森序列序列。

②.由递归辛普森序列序列得到递归布尔公式序列。

③.通过理查森改进提高误差项的阶数。

4.用自适应积分方法计算
图4.自适应积分算法流程图
在辛普森公式基础上，将区间再进行划分，即为自适应积分。

3. 实验结果及分析
真实值S = 0.19971466216144
1. 用组合梯形公式M=10计算。

面积近似值S1 =0.16965032127666。

误差error1=0.03006434088479。

2.用组合辛普森公式M=5计算
4.660686426147481*10-。

面积近似值S2 =0.19966805529718。

误差error2=5 3.用龙贝格积分计算。

表1.龙贝格积分表
面积近似值S3 =0.19971378242654。

误差error3=5
2.254986251026825*10-。

4.用自适应积分方法计算
面积近似值S4 =0.19971391278871。

误差error4=6
3.238553837777504*10-。

4. 结论
自适应计分方法和龙贝格积分法都可以得到较高的精确度。

组合梯形和组合辛普森可以增大M值，得到较高的精确度。

附件（代码）
1.
f=inline('sin(4*x)*exp(-2*x)','x');%定义函数
A=0;t=0;S=0;
for k=1:9
t=k*0.3;
A=A+0.3*feval(f,t);
end
S=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;
end
for 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))); end
RS=int(sin(4*x)*exp(-2*x),0,3);%真实值
S
error=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 is the quadrature value
% -err is the error estimate
% -h is the smallest step size uesd
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
while((err>tol)&(J<n))|(J<4)
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)+(R(J+1,K)-R(J,K))/(4^K-1);
end
err=abs(R(J,J)-R(J+1,K+1));
end
quad=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 tol1]
h=(b0-a0)/2;
C=zeros(1,3);
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 is the quadrature value
% -err is the error estimate
% -Initialize values
SRmat=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)
%Bisect interval,apply Simpson'rule
%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 (err<tol)
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
end
end
quad=sum(SRmat(:,4));
err=sum(abs(SRmat(:,5)));
SRmat=SRmat(1:m,1:6);。