matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

9.3.6 牛顿-科茨(Newton-Cotes)数值积分和误差分析的MATLAB程序

9.3.6  牛顿-科茨(Newton-Cotes)数值积分和误差分析的MATLAB程序

 

(一)  估计误差的MATLAB程序

计算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

 

例9.3.14  用求截断误差公式的MATLAB主程序,求计算定积分d的近似值的阶牛顿–科茨公式的截断误差公式.

解 输入求阶牛顿–科茨公式的截断误差公式的程序

>> 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 =          RNC1 =

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

n =          RNC2 =

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

n =          RNC3 =

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

n =          RNC4 =

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

n =          RNC8 =

     8          35065906189543/6926923254988800*(1/8*b-1/8*a)^11*fxn2

 

例9.3.15  用MATLAB软件估计用5、6阶牛顿–科茨公式计算定积分ed的截断误差公式和误差限,取15位小数.

解  输入求阶牛顿–科茨公式的截断误差公式和被积函数的6,8阶导函数的程序

>> n=5;RNC5=ncE(n), n=6;RNC6=ncE(n)

syms x

y=exp(sin(x)); yx6=diff(y,x,6), yx8=diff(y,x,8)

运行后输出被积函数的6,8阶导函数(略)和阶牛顿–科茨公式的截断误差公式如下

          RNC5 =                               RNC6 =

275/12096*(1/5*b-1/5*a)^7*fxn1      9/1400*(1/6*b-1/6*a)^9*fxn2

然后再输入误差估计程序

>>a=0;b=pi/2; h=pi/40; x=0:0.00001:pi/2;

yx6=-sin(x).*exp(sin(x))+16*cos(x).^2.*exp(sin(x))-15*sin(x).^2.*exp(sin(x))+75*sin(x).*cos(x).^2.*exp(sin(x))-20*cos(x).^4.*exp(sin(x))-15*sin(x).^3.*exp(sin(x))+45*sin(x).^2.*cos(x).^2.*exp(sin(x))-15*sin(x).*cos(x).^4.*exp(sin(x))+cos(x).^6.*exp(sin(x));

yx8=cos(x).^8.*exp(sin(x))-756*sin(x).*cos(x).^2.*exp(sin(x))-1260*sin(x).^2.*cos(x).^2.*exp(sin(x))+630*sin(x).*cos(x).^4.*exp(sin(x))-420*sin(x).^3.*cos(x).^2.*exp(sin(x))+210*sin(x).^2.*cos(x).^4.*exp(sin(x))-28*sin(x).*cos(x).^6.*exp(sin(x))-56*cos(x).^6.*exp(sin(x))+sin(x).*exp(sin(x))+63*sin(x).^2.*exp(sin(x))+210*sin(x).^3.*exp(sin(x))+105*sin(x).^4.*exp(sin(x))-64*cos(x).^2.*exp(sin(x))+336*cos(x).^4.*exp(sin(x));

myx6=max(yx6); myx8=max(yx8); RNC5 =275/12096*(1/5*b-1/5*a)^7*myx6

RNC6 =9/1400*(1/6*b-1/6*a)^9*myx8

运行后屏幕显示误差限如下

  RNC5 =                           RNC6 =

    3.625385713339320e-004               3.826183594914823e-005

 

(二) 计算科茨系数和求截断误差公式的MATLAB程序

计算n阶科茨系数和求截断误差公式的MATLAB主程序

function [Cn, RNCn]=newcotE(n)

syms t a b M,Fz=zeros(1,n+1);

Cn=zeros(1,n+1); su=t;k=1;m=1;m0=1;

  for u=1:n

su1=su*(t-u);m01=m0*u; su=su1;m0=m01;

end

su;m0; f1=su/(t-0); intf1=int(f1,t,0,n);

y= double(intf1);

Cn(1)=((-1)^(n-0)*y)/(n*m0); k=1;m=1;

for j=1:n

k1=k*j; m1=m*(n-j); f=su/(t-j);

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

Cn(j+1)=((-1)^(n-j)*y)/(n*k1*m1);

warning off MATLAB:divideByZero

end

 fn=su/(t-n); intfn=int(fn,t,0,n);

y= double(intfn);Cn(n+1)=y/(n*m0);

Cn; 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);

RNCn=(((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);

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

end

 

例9.3.16  用计算n阶科茨系数和求截断误差公式的MATLAB主程序,计算定积分I=d阶牛顿–科茨公式的系数和截断误差公式.

解 (1)先求阶牛顿–科茨公式的系数和截断误差公式.输入程序

>>  n1=1,[Cn1,RNCn1]=newcotE(n1)

n2=2,[Cn2,RNCn2]=newcotE(n2)

n3=3,[Cn3,RNCn3]=newcotE(n3)

运行后屏幕显示阶牛顿–科茨公式的系数Cn1, Cn2, Cn3和截断误差公式RNCn1, RNCn2, RNCn3如下

n1 =

      1     

Cn1 =

     1/2          1/2    

RNCn1 =

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

n2 =

      2     

Cn2 =

     1/6          2/3          1/6    

RNCn2 =

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

n3 =

      3     

Cn3 =

     1/8          3/8          3/8          1/8    

RNCn3 =

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

 

(三) 计算牛顿–科茨公式的MATLAB程序

 

调用格式一:quad8(‘fun’,a,b)

调用格式二:quad8(‘fun’,a,b,tol)

    调用格式三: quad8(‘fun’,a,b, tol,TRACE)

调用格式四:quad8(‘fun’,a,b, tol,TRACE,P1,P2,...)

 

例9.3.17  用梯形公式、辛普森和牛顿–科茨求积公式计算定积分ed,取精度,作出它们的积分图,并与精确值进行比较;

 

解  (1)用梯形求积公式计算定积分. 输入程序

>> h=pi/500; x=0:h:pi/2; y=exp(sin(x));

zt=trapz(x,y), ztc=cumtrapz(x,y), plot(x, ztc,'ro')

运行后屏幕显示用函数trapz 和cumtrapz分别计算结果zt、ztc分别如下

zt =

   3.10437572798742

ztc =

  Columns 1 through 3

                  0   0.00630298652792   0.01264569951380

………………………………………………………………………..

Columns 250 through 251

   3.08729642810745   3.10437572798742

(2)用辛普森求积公式计算定积分. 输入程序

>> syms x

L= inline(' exp(sin(x))');

[QS,FCNTS] =quad(L,0, pi/2,1.e-4,2)

运行后屏幕显示用辛普森求积公式计算定积分的值QS和递归次数FCNTS分别如下

QS =                           FCNTS =

   3.10438133817254                    13

(3)用牛顿–科茨求积公式计算定积分. 在MATLAB6.5中输入程序

>> syms x

L= inline(' exp(sin(x))');

[Q8,FCNT8] = quad8(L,0, pi/2,1.e-4,3)

运行后屏幕显示用牛顿–科茨求积公式计算定积分的值Q8和递归次数FCNTS分别如下

Q8 =                              FCNT8 =

        3.10437901785555                      33

(4)输入求定积分的精确值的程序

>> syms x

y=exp(sin(x)); F=int(y,0,pi/2), Fj=double(F)

wzt=abs( Fj- zt), wQS = abs(Fj- QS), wQ8 = abs(Fj- Q8)

运行后屏幕显示计算的定积分值F和近似值Fj,梯形公式、辛普森和牛顿–科茨求积公式计算定积分的值与Fj的绝对误差wztc, wQS和wQ8如下

Warning: Explicit integral could not be found.

> In C:\MATLAB6p5p1\toolbox\symbolic\@sym\int.m at line 58

F =

int(exp(sin(x)),x = 0 .. 1/2*pi)

Fj =                             wzt =

  3.10437901785556                    3.289868133471430e-006

wQS =                           wQ8 =

  2.320316987436399e-006             7.993605777301127e-015

(5)输入画图程序

>> syms x

F=int(exp(sin(x)),0,pi/2),Fj=double(F),plot(Fj,'ro'),

hold on

L= inline(' exp(sin(x))');

Q=quad(L,0, pi/2,1.e-4,2),plot(Q,'g*')

hold off,hold on, h=pi/40;

x=0:h:pi/2; y=exp(sin(x));

ztc=cumtrapz(x,y),

plot(x, ztc,'mp'), hold off

hold on, [Q8,FCNT8] = quad8(L,0, pi/2,1.e-4,3), hold off

grid, xlabel('自变量 X'), ylabel('因变量 Y')

title('牛顿–科茨公式和梯形公式计算数值积分')

legend('精确值', ' 辛普森公式计算定积分','梯形公式计算定积分','牛顿–科茨公式计算定积分')

运行后屏幕显示图形(略).

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3090306.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1