1 梯形数值积分的MATLAB主程序

function T=rctrap(fun,a,b,m)

%fun 函数，a 积分上限 b积分下限 m 递归次数

n=1;h=b-a; T=zeros(1,m+1); x=a;

T(1)=h*(feval_r(fun,a)+feval_r(fun,b))/2;

for i=1:m

h=h/2; n=2*n; s=0;

for k=1:n/2

x=a+h*(2*k-1); s=s+feval_r(fun,x);

end

T(i+1)=T(i)/2+h*s;

end

T=T(1:m);

e.g

运行后屏幕显示 精确值Fs，用rctrap计算的递归值T和T与精确值Fs的绝对误差wT

>> fun=@(x)

exp((-x^.2./2)./(sqrt(2*pi)))

T=rctrap(fun,0,pi/2,14), syms t

fi=int(exp((-t^2)/2)/(sqrt(2*pi)),t,0, pi/2);

Fs= double(fi), wT= double(abs(fi-T))

fun =

@(x)exp((-x^.2./2)./(sqrt(2*pi)))

T =

Columns 1 through 7

1.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109

Columns 8 through 14

1.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101

Fs =

0.4419

wT =

Columns 1 through 7

0.9749 0.9159 0.8894 0.8776 0.8723 0.8700 0.8690

Columns 8 through 14

0.8686 0.8684 0.8683 0.8683 0.8683 0.8682 0.8682

>>

2 复合辛普森(Simpson)数值积分的MATLAB主程序

function y=comsimpson(fun,a,b,n)

% fun 函数 a 积分上限 b积分下限 n 分割小区间数

z1=feval_r(fun,a)+ feval_r(fun,b);m=n/2;

h=(b-a)/(2*m); x=a;

z2=0; z3=0; x2=0; x3=0;

for k=2:2:2*m

x2=x+k*h; z2= z2+2*feval_r(fun,x2);

end

for k=3:2:2*m

x3=x+k*h; z3= z3+4*feval_r(fun,x3);

end

y=(z1+z2+z3)*h/3;

由于Matlab自带了 quad就是这个算法 所以比较少自己编

3 龙贝格数值积分的MATLAB主程序

function [RT,R,wugu,h]=romberg(fun,a,b,

wucha,m)

%fun被积函数 a，b积分上下限 wucha两次相邻迭代绝对差值 m 龙贝格积分表最大行数

%RT 龙贝格积分表 R 数值积分结果 wucha 误差估计 h 最小步长

n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4);

RT(1,1)=h*(feval_r(fun,a)+feval_r(fun,b))/2;

while((wugu>wucha)&(k

k=k+1; h=h/2; s=0;

for j=1:n

x=a+h*(2*j-1); s=s+feval_r(fun,x);

end

RT(k+1,1)= RT(k,1)/2+h*s; n=2*n;

for i=1:k

RT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);

end

wugu=abs(RT(k+1,k)-RT(k+1,k+1));

end

R=RT(k+1,k+1);

e.g.

>> F=inline('1./(1+x)');

[RT,R,wugu,h]=romberg(F,0,1.5,1.e-8,13)

syms x

fi=int(1/(1+x),x,0,1.5); Fs=double(fi),

wR=double(abs(fi-R)), wR1= wR - wugu

RT =

1.0500 0 0 0 0 0

0.9536 0.9214 0 0 0 0

0.9260 0.9168 0.9165 0 0 0

0.9187 0.9163 0.9163 0.9163 0 0

0.9169 0.9163 0.9163 0.9163 0.9163 0

0.9164 0.9163 0.9163 0.9163 0.9163 0.9163

R =

0.9163

wugu =

2.9436e-011

h =

0.0469

Fs =

0.9163

wR =

9.8007e-011

wR1 =

6.8571e-011

>>

4 复合梯形法function [I,step]

= CombineTraprl(f,a,b,eps)

%f 被积函数

%a,b 积分上下限

%eps 精度

%I 积分结果

%step 积分的子区间数

if(nargin ==3)

eps=1.0e-4;

end

n=1;

h=(b-a)/2;

I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;

while abs(I2-I1)>eps

n=n+1;

h=(b-a)/n;

I1=I2;

I2=0;

for

i=0:n-1

x=a+h*i;

x1=x+h;

I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));

end

end

I=I2;

step=n;

5 辛普森法

function [I,step] = IntSimpson(f,a,b,type,eps)

%type = 1 辛普森公式

%type = 2 辛普森3/8公式

%type = 3 复合辛普森公式

if(type==3 && nargin==4)

eps=1.0e-4; %缺省精度为0.0001

end

I=0;

switch type

case

1,

I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+...

4*subs(sym(f),findsym(sym(f)),(a+b)/2)+...

subs(sym(f),findsym(sym(f)),b));

step=1;

case

2,

I=((b-a)/8)*(subs(sym(f),findsym(sym(f)),a)+...

3*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+ ...

3*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+subs(sym(f),findsym(sym(f)),b));

step=1;

case

3,

n=2;

h=(b-a)/2;

I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;

while abs(I2-I1)>eps

n=n+1;

h=(b-a)/n;

I1=I2;

I2=0;

for i=0:n-1

x=a+h*i;

x1=x+h;

I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+...