2024年5月26日发(作者：浦代卉)

1．用求截断误差公式的MATLAB主程序，求计

算定积分



b

a

f(x)

d

x

的近似值的

1,2,3,4,8

阶牛顿–

科茨公式的截断误差公式。

Matlab程序如下

n=1, RNC1=ncE(n), n=2, RNC2=ncE(n),

n=3, RNC3=ncE(n)

n=4, RNC4=ncE(n), n=8, RNC8=ncE(n)

运行后屏幕显示结果如下

n= 1 RNC1 =1/12*(b-a)^3*fxn1

n=2 RNC2 =1/90*(1/2*b-1/2*a)^5*fxn2

n = 3 RNC3=3/80*(1/3*b-1/3*a)^5*fxn1

n =4 RNC4=8/945*(1/4*b-1/4*a)^7*fxn2

n=8 RNC8=

35/6926923254988800*(1/8*b-1/8*a)

^11*fxn2

其中计算n阶牛顿-科茨的公式的截断误差

公式的MATLAB主程序

function RNC=ncE(n)

suk=1; p=n/2-fix(n/2);

if p==0

for k=1:n+2

suk=suk*k;

end

suk; syms t a b fxn2,su=t^2;

for u=1:n

su=su*(t-u);

end

su; intf=int(su,t,0,n); y=double(intf);

RNC= (((b-a)/n)^(n+3))*fxn2*abs(y)/

suk;

else

for k=1:n+1

suk=suk*k;

end

suk; syms t a b fxn1,su=t;

for u=1:n

su=su*(t-u);

end

su; intf=int(su,t,0,n);

y=double(intf);

RNC=

(((b-a)/n)^(n+2))*fxn1*abs(y)/ suk;

end

2．分别利用复化梯形公式和复化Simpson公

式计算定积分

I



2

0

1e

x

dx

取

n2,4,8,16

，精确值为

I4.006994

。

复化梯形公式MATLAB程序如下

clear;

Iexact=4.006994;

a=0;b=2;

fprintf('n n I Errorn');

n=1;

for k=1:4

n=2*n;

h=(b-a)/n;i=1:n+1;

x=a+(i-1)*h;f=sqrt(1+exp(x));

I=trapez_v(f,h);

fprintf('%3.0f %10.5f %10.5fn',n,I,Iexact-I)

;

end

其中复化梯形求积法的M文件如下

trapez_v.m

function I=trapez_v(f,h)

I=h*[sum(f)-(f(1)+f(length(f)))/2];

结果：n =2 I =4.08358 误差为-0.07659

n =4 I =4.02619 误差为

-0.01919

n =8 I =4.01180 误差为

-0.00480

n =16 I =4.00819 误差为

-0.00120

复化Simpson求积公式MATLAB程序如下：

clear;

Iexact=4.006994;

a=0;b=2;

fprintf('n n I Errorn');

n=1;

for k=1:4,n=2*n

h=(b-a)/n;i=1:n+1;

x=a+(i-1)*h;f=sqrt(1+exp(x));

I=(h/3)*(f(1)+4*sum(f(2:2:n))+f(n+1));

if n>2

I=I+(h/3)*2*sum(f(3:2:n));

end

fprintf('%3.0f %10.5f %10.5fn',n,I,Iexact-I);

end

结果n =2 I=4.00791 误差为-0.00092

n =4 I=4.00705 误差为-0.00006

n =8 I=4.00700 误差为-0.00000

n =16 I=4.00699 误差为-0.00000

3．计算无穷限广义积分

I

1



x

2







edx

，其精确值

I1

。

用-10,10替换积分限，则有

I

1



10

10

e

x

2



dx

利用复化梯形公式，调用trapez_v.m文件，取

n=10,20的MATLAB程序如下：

I=1;

a=-10;b=10;

fprintf('n h n I1n');

n=0;

for k=1:2;n=n+20;

h=(b-a)/n;i=1:n+1;

x=a+(i-1)*h;f=exp(-x.^2);

I=trapez_v(f,h);

I1=1/sqrt(pi)*I;

fprintf('%3.1f %10.0f %10.8fn',h,n,I1);

end

结果 h=1.0 n=20 I1=1.00010345

h=0.5 n=40 I1=1.00000000

4．取精度为



10

8

，分别用

T

j1,j1

T

j1,j





和

T

j1,j1

T

j,j





作为计算停止的条件，用Romberg程序计算

I



1.5

1

0

1x

d

x

，取精度为

10

-8

,并与精确值

比较。然后取精度为

10

-7

，用

T

j1,j1

T

j1,j





作为计算停止的条件，观

察用龙贝格求积公式计算的结果与精确值的绝对

误差是否满足

10

-7

。

MATLAB程序如下

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

其中龙贝格积分的MATLAB程序

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

wucha,m)

n=1;h=b-a; wugu=1; x=a;k=0;

RT=zeros(4,4);

RT(1,1)=h*(feval(fun,a)+feval(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(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);

结果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

2024年5月26日发(作者：浦代卉)

1．用求截断误差公式的MATLAB主程序，求计

算定积分



b

a

f(x)

d

x

的近似值的

1,2,3,4,8

阶牛顿–

科茨公式的截断误差公式。

Matlab程序如下

n=1, RNC1=ncE(n), n=2, RNC2=ncE(n),

n=3, RNC3=ncE(n)

n=4, RNC4=ncE(n), n=8, RNC8=ncE(n)

运行后屏幕显示结果如下

n= 1 RNC1 =1/12*(b-a)^3*fxn1

n=2 RNC2 =1/90*(1/2*b-1/2*a)^5*fxn2

n = 3 RNC3=3/80*(1/3*b-1/3*a)^5*fxn1

n =4 RNC4=8/945*(1/4*b-1/4*a)^7*fxn2

n=8 RNC8=

35/6926923254988800*(1/8*b-1/8*a)

^11*fxn2

其中计算n阶牛顿-科茨的公式的截断误差

公式的MATLAB主程序

function RNC=ncE(n)

suk=1; p=n/2-fix(n/2);

if p==0

for k=1:n+2

suk=suk*k;

end

suk; syms t a b fxn2,su=t^2;

for u=1:n

su=su*(t-u);

end

su; intf=int(su,t,0,n); y=double(intf);

RNC= (((b-a)/n)^(n+3))*fxn2*abs(y)/

suk;

else

for k=1:n+1

suk=suk*k;

end

suk; syms t a b fxn1,su=t;

for u=1:n

su=su*(t-u);

end

su; intf=int(su,t,0,n);

y=double(intf);

RNC=

(((b-a)/n)^(n+2))*fxn1*abs(y)/ suk;

end

2．分别利用复化梯形公式和复化Simpson公

式计算定积分

I



2

0

1e

x

dx

取

n2,4,8,16

，精确值为

I4.006994

。

复化梯形公式MATLAB程序如下

clear;

Iexact=4.006994;

a=0;b=2;

fprintf('n n I Errorn');

n=1;

for k=1:4

n=2*n;

h=(b-a)/n;i=1:n+1;

x=a+(i-1)*h;f=sqrt(1+exp(x));

I=trapez_v(f,h);

fprintf('%3.0f %10.5f %10.5fn',n,I,Iexact-I)

;

end

其中复化梯形求积法的M文件如下

trapez_v.m

function I=trapez_v(f,h)

I=h*[sum(f)-(f(1)+f(length(f)))/2];

结果：n =2 I =4.08358 误差为-0.07659

n =4 I =4.02619 误差为

-0.01919

n =8 I =4.01180 误差为

-0.00480

n =16 I =4.00819 误差为

-0.00120

复化Simpson求积公式MATLAB程序如下：

clear;

Iexact=4.006994;

a=0;b=2;

fprintf('n n I Errorn');

n=1;

for k=1:4,n=2*n

h=(b-a)/n;i=1:n+1;

x=a+(i-1)*h;f=sqrt(1+exp(x));

I=(h/3)*(f(1)+4*sum(f(2:2:n))+f(n+1));

if n>2

I=I+(h/3)*2*sum(f(3:2:n));

end

fprintf('%3.0f %10.5f %10.5fn',n,I,Iexact-I);

end

结果n =2 I=4.00791 误差为-0.00092

n =4 I=4.00705 误差为-0.00006

n =8 I=4.00700 误差为-0.00000

n =16 I=4.00699 误差为-0.00000

3．计算无穷限广义积分

I

1



x

2







edx

，其精确值

I1

。

用-10,10替换积分限，则有

I

1



10

10

e

x

2



dx

利用复化梯形公式，调用trapez_v.m文件，取

n=10,20的MATLAB程序如下：

I=1;

a=-10;b=10;

fprintf('n h n I1n');

n=0;

for k=1:2;n=n+20;

h=(b-a)/n;i=1:n+1;

x=a+(i-1)*h;f=exp(-x.^2);

I=trapez_v(f,h);

I1=1/sqrt(pi)*I;

fprintf('%3.1f %10.0f %10.8fn',h,n,I1);

end

结果 h=1.0 n=20 I1=1.00010345

h=0.5 n=40 I1=1.00000000

4．取精度为



10

8

，分别用

T

j1,j1

T

j1,j





和

T

j1,j1

T

j,j





作为计算停止的条件，用Romberg程序计算

I



1.5

1

0

1x

d

x

，取精度为

10

-8

,并与精确值

比较。然后取精度为

10

-7

，用

T

j1,j1

T

j1,j





作为计算停止的条件，观

察用龙贝格求积公式计算的结果与精确值的绝对

误差是否满足

10

-7

。

MATLAB程序如下

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

其中龙贝格积分的MATLAB程序

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

wucha,m)

n=1;h=b-a; wugu=1; x=a;k=0;

RT=zeros(4,4);

RT(1,1)=h*(feval(fun,a)+feval(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(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);

结果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