积分方程的数值计算方法

S=0;X=zeros(1,11);h=0.3;
forj=1:11
X(j)=(j-1)*h;
end
fork=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);%真实值
A=0;t=0;S=0;
fork=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
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
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;
% -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
2.用组合辛普森公式M=5计算
面积近似值S2=0.19966805529718。误差error2= 。
3.用龙贝格积分计算。
表1.龙贝格积分表
J
梯形公式
辛普森公式
布尔公式
0
-0.00199504693265
0
0
0
0
0
1
-0.02186444123413
-0.02848757266796
0
0
0
0
2
0.01611754400127
0.19974037084997
0.19969123256403
0
5
0.19679137798386
0.19972090351497
0.19971665225786
0.19971386435224
0.19971376040519
0.19971378242654
面积近似值S3=0.19971378242654。误差error3= 。
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;
forK=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);
2.用组合辛普生公式M=5计算。
图2. 组合辛普森算法流程图
将积分区间 划分为宽度为 的2M个子区间，再由辛普森公式将各区间的面积求和即可得到近似面积。
3.用龙贝格积分计算。
图3.龙贝格积分算法流程图
.由递归梯形公式序列得到递归辛普森序列序列。
.由递归辛普森序列序列得到递归布尔公式序列。
.通过理查森改进提高误差项的阶数。
% -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;
forp=1:M
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;
0.02877820574640
0.03259592430736
0
0
0
3
0.15265718844742
0.19817040326280
0.20946321643056
0.21227063376585
0
0
4
0.18800280139054
0.19978467237158
0.19989229031217
end
end
end
end
quad=sum(SRmat(:,4));
err=sum(abs(SRmat(:,5)));
SRmat=SRmat(1:m,1:6);
实验七 积分方程的数值计算方法
1. 实验描述
计算 定积分的近似值，起始容差
1.用组合梯形公式M=10计算。
2.用组合辛普生公式M=5计算。
3.用龙贝格积分计算。
4.用自适应积分方法计算。
2.实验内容
1. 用组合梯形公式M=10计算。
图1. 组合梯形算法流程图
将积分区间 划分为宽度为 的M个子区间，再将各区间的面积求和即可得到近似面积。
SRvec=srule(f,a,b,tol);
SRmat(1,1:6)=SRvec;
m=1;
state=iterating;
while(state==iterating)
n=m;
forj=n:-1:1
p=j;
SR0vec=SRmat(p,:);
err=SR0vec(5);
tol=SR0vec(6);
4.用自适应积分方法计算
面积近似值S4=0.19971391278871。误差error4= 。
4.结论
自适应计分方法和龙贝格积分法都可以得到较高的精确度。
组合梯形和组合辛普森可以增大M值，得到较高的精确度。
附件（代码）
1.
f=inline('sin(4*x)*exp(-2*x)','x');%定义函数
4.用自适应积分方法计算
图4.自适应积分算法流程图
在辛普森公式基础上，将区间再进行划分，即为自适应积分。
3.实验结果及分析
真实值S= 0.19971466216144
1. 用组合梯形公式M=10ຫໍສະໝຸດ Baidu算。
面积近似值S1=0.16965032127666。误差error1=0.03006434088479。
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)
%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);
4.
functionZ=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