积分方程的数值计算技巧

实验七积分方程的数值计算方法

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

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;