




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、第九章 数值积分积分是非常重要的数学工具(微分方程、概率论等的基础),另一方面它们在实际问题中也有着许多直接的应用.函数的积分计算可分为数值积分和符号积分两类方法.数值积分是求积分的近似值的一类近似计算的方法.对于定积分d,如果对应的不定积分d不易求出,或者根本不能表示为初等函数时,那么就只能用数值积分的方法求其近似值了.例如,因为不定积分ed,d,ed,d无法由初等函数表示,所以计算这种类型的定积分只能用数值方法.至于由离散数据或图形表示的函数的定积分,理所当然地属于数值积分的范畴.符号积分是指求积分的精确值.本章重点介绍数值积分及其用MATLAB计算的方法.为了将数值积分的值与精确值比较,
2、在第一节中简单介绍用MATLAB 作符号积分的方法.9.1 积分的MATLAB符号计算计算机软件MATLAB的系统提供了符号计算定积分的函数int.下面主要介绍用MATLAB软件进行符号计算定积分、变上(下)限积分及其应用.9.1.1 定积分的MATLAB符号计算如果函数在区间上连续,且是的一个原函数,则函数在区间上的定积分为d.我们可以用MATLAB软件中的函数int对定积分进行符号计算,它的调用格式和功能如表91所示.表91 函数int求定积分的调用格式和功能函数int求定积分的调用格式功 能F=int(S,a,b)输入量:S可以是被积函数,也可以是由几个被积函数,构成的矩阵,积分上、下限
3、a 和 b可以是数,也可以是符号标量,积分变量是.输出量:F是S关于积分变量的从a 到 b的定积分.F=int(S,t,a,b)输入量:S、a 和 b同上,是积分变量.输出量:F是S关于积分变量从a 到 b的定积分.例9.1.5 由,所围成的平面区域D.求平面区域D的面积S.解 (1) 画出由,所围成的封闭平面区域的图形,输入作函数图形的程序>> x=-1:0.001:2; F1= sin(x); F2=cos(x);plot(x ,F1,'b-',x ,F2,'g-'), axis(-1,pi/4+1,-1.3,1.3),xlabel('x
4、'), ylabel('y'),title('y=sinx , y=cosx 和x=-0.5及x=1.5所围成的平面区域的图形')运行后屏幕显示图91.图91 由,所围成的平面区域图形(2)求平面区域D的面积S.解方程组得在内的两条曲线交点的横坐标,所求面积为dd.输入计算面积S的程序 >> syms x f1= cos(x)-sin(x); f2=-f1; S1 =int(f1,x,-0.5,pi/4); S2=int(f2,x, pi/4,1.5); S=S1+S2,Sj= double (S)运行后屏幕显示计算面积的值 S 及其近似值
5、Sj 如下S =2*2(1/2)+sin(1/2)-cos(1/2)-sin(3/2)-cos(3/2)Sj = 1.362037913188269.1.2 变限积分的MATLAB符号计算例9.1.7 已知ed,求.解 因为edteded所以输入程序 >> syms x tF1=int(exp(t)*sin(2+sqrt(t3),x,0);F2=int(exp(t)*sin(2+sqrt(t3),0,x2);Fi= F1+ F2;dF=diff(Fi)运行后屏幕显示计算变限积分的导数如下dF =-exp(x)*sin(2+(x3)(1/2)+2*x*exp(x2)*sin(2+(x
6、6)(1/2)9.2 数值积分的思想及其MATLAB程序9.2.3 矩形公式的MATLAB程序用矩形公式计算数值积分可用下面函数sum和cumsum实现,调用格式如下:(一) 函数sum的调用格式函数sum的两种调用格式如下:调用格式一:sum(X)如果输入向量X,则输出向量X的所有元素的和,可用于按矩形公式(9.3)和(9.4)计算积分;如果输入矩阵X,则输出矩阵X的每列向量的和.调用格式二:sum (X,DIM)输入N-D 数组X,输出为X的每个DIM向量的和.例9.2.2 用MATLAB和矩形公式(9.3)、(9.4)计算ed,并与精确值比较.解 将分成20等分,步长为,输入程序>
7、> h=pi/40; x=0:h:pi/2; y=exp(sin(x);z1=sum(y(1:20)*h, z2=sum(y(2:21)*h,运行后屏幕显示矩形公式计算结果分别如下 z1 = z2 = 3.0364 3.1713即矩形公式(9.3)计算结果3.036 4, 矩形公式(9.4)计算结3.171 3.求定积分的精确值,输入程序>> syms x F=int(exp(sin(x),x,0, pi/2), Fs= double (F), wz1=abs( Fs-z1), wz2= abs( Fs-z2)运行后屏幕显示定积分的精确值Fs和与用矩形公式(9.3)、(9.4
8、)计算结果的绝对误差wz1、wz2分别如下Warning: Explicit integral could not be found.> In C:MATLAB6p5p1toolboxsymbolicsymint.m at line 58F = Fs =int(exp(sin(x),x = 0 . 1/2*pi) 3.1044wz1 = wz2 = 0.0680 0.0670由此可见,该定积分用矩形公式(9.4)计算结果比矩形公式(9.3)的绝对误差小.(二) 函数cumsum的调用格式函数cumsum的调用格式如下:调用格式一:cumsum(X) 输入向量X,输出向量X的元素依次累加和
9、,可用于按矩形公式(9.3)、(9.4)计算积分; 输入矩阵X,输出矩阵为X的每列向量的依次累加和.调用格式二:cumsum (X,DIM)输入N-D 数组X,输出为X的每个DIM向量的依次累加和.例9.2.4 用MATLAB的函数sum 和 cumsum及矩形公式(9.3)、(9.4)计算ed,并与精确值比较.解 将分成20等分,步长为,输入程序如下(注意sum 和 cumsum的用法)>> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x); z1=sum(y(1:20)*h,z2=sum(y(2:21)*h, z=cumsum(y); z11=z(2
10、0)*h, z12=(z(21)-z(1)*h,运行后屏幕显示计算结果分别如下z1 = z2 = z11 = z12 = 0.3873 0.4036 0.3873 0.4036求定积分的精确值,输入程序>> syms x F=int(exp(-x)*sin(x),x,0, pi/2)Fs= double (F) ,wz1=abs( Fs-z1), wz2= abs( Fs-z2)运行后屏幕显示定积分的精确值Fs和用矩形公式(9.3)、(9.4)计算结果的绝对误差wz1、wz2分别如下F = Fs =1/2*(-1+exp(pi)(1/2)/exp(pi)(1/2) 0.3961wz
11、1 = wz2 = 0.0088 0.0075由此可见,用矩形公式(9.4)计算该定积分的结果比用矩形公式(9.3)的结果的绝对误差小.9.3 插值型数值积分及其MATLAB 程序如果被积函数和积分变量是离散的数据X,Y或数表,则常用的处理方法是利用各种插值函数(例如,样条插值函数、牛顿插值函数、拉格朗日插值多项式等)构造数值积分公式.本节主要介绍利用拉格朗日多项式和样条插值函数等构造数值积分公式的方法、梯形公式、辛普森(Simpson)公式和一般的Newton-Cotes公式及其误差分析,三者之间的关系,并通过实例说明用MATLAB软件计算数值积分的方法.9.3.2 梯形公式的MATLAB程
12、序(一) 根据梯形公式和估计误差公式自己编写MATLAB程序计算定积分例9.3.2 分别取,用梯形公式计算定积分ed,并与精确值比较.然后观察对计算结果的有效数字和绝对误差的影响.解 根据(9.5a)式编写并输入如下程序>> h=pi/8000; a=0;b=pi/2; x=a:h:b;n=length(x), y=exp(sin(x);z1=(y(1)+y(n)*h/2; z2=sum(y(2:n-1)*h; z8000=z1+z2,syms tf=exp(sin(t); intf=int(f,t,a,b), Fs=double(intf),Juewucha8000=abs(z8
13、000-Fs)运行后屏幕显示取时,积分区间上等距节点的个数,用梯形公式计算定积分的值z8000和精确值intf的近似值Fs及其绝对误差Juewucha8000如下n = z8000 = 4001 3.10437900500451Warning: Explicit integral could not be found.> In C:MATLAB6p5toolboxsymbolicsymint.m at line 58intf = Fs =int(exp(sin(t),t = 0 . 1/2*pi) 3.10437901785556Juewucha8000 = 1.285104200832
14、166e-008然后再分别取,运行后屏幕依次显示用梯形公式计算定积分的值z800、z80和与精确值intf的近似值Fs的绝对误差Juewucha 800、Juewucha 80如下n800 = z800 = Juewucha800 = 401 3.10437773275082 1.285104739068288e-006n80 = z80= Juewucha80 =41 3.10425050738255 1.285104730022191e-004 (二) 用函数trapz计算定积分用MATLAB软件和梯形公式计算数值积分和估计误差限除了用例9.3.2和例9.3.3自己编写的MATLAB计算程
15、序外,MATLAB系统提供了两种计算程序:一种是梯形数值积分(Trapezoidal doubleal integration.)的程序trapz.m;另一种是累加梯形数值积分(Cumulative trapezoidal doubleal integration.)的程序cumtrapz.m,它们的调用格式如下:调用格式一:Z =trapz(Y)输入Y,输出为按梯形公式(9.5a)计算的Y的积分的近似值(单位步长).(1)如果输入Y是向量,则Z是Y的积分;(2)如果输入Y是矩阵,则行向量Z是Y的每列的积分.调用格式二:Z =trapz(X,Y)输入X,Y为同长度的数组,输出为按梯形公式(9.
16、5a)计算Y对X的积分(步长不一定相等).调用格式三:Z = trapz (X,Y,DIM) 或 trapz (Y,DIM)输出为按梯形公式(9.5a)计算Y 的每个DIM的积分,其中X的长度必须和size(Y,DIM)相同.(三) 用函数cumtrapz计算定积分调用格式一:Z =cumtrapz (Y)输入Y,输出用梯形方法计算的Y的累加积分(Cumulative integral)的近似值(单位步长).(1)如果输入Y是向量,则Z是Y的累加积分;(2)如果输入Y是矩阵,则行向量Z是Y的每列的累加积分.调用格式二:Z =cumtrapz (X,Y)输入X,Y为同长度的数组,输出用梯形方法计
17、算的Y对X的累加积分(步长不一定相等).调用格式三:Z = cumtrapz (X,Y,DIM) 或 cumtrapz (Y,DIM)输出用梯形方法计算的Y 的每个DIM的累加积分,其中X的长度必须和size(Y,DIM)相同. 例9.3.4 用MATLAB的函数trapz 和cumtrapz分别计算ed,精确到,并与矩形公式(9.3),(9.4)比较.解 将分成20等分,步长为,输入程序如下(注意trapz(y)是单位步长, trapz(y)*h=trapz(x,y)):>> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x);z1=sum(y(1:20
18、)*h, z2=sum(y(2:21)*h, z=(z1+z2)/2z3=trapz(y)*h, z3h=trapz(x,y), z3c=cumtrapz(y)*h,运行后屏幕显示用矩形公式(9.3),(9.4)计算结果z1、z2和二者的平均数z、函数trapz 和cumtrapz分别计算结果z3、z3c分别如下z1 = z2 = z = z3 = z3h =0.3873 0.4036 0.3954 0.3954 0.3954z3c = Columns 1 through 7 0 0.0028 0.0109 0.0234 0.0395 0.0586 0.0798 Columns 8 throu
19、gh 14 0.1028 0.1270 0.1519 0.1771 0.2023 0.2273 0.2517 Columns 15 through 21 0.2755 0.2983 0.3201 0.3408 0.3602 0.3785 0.3954显然,函数trapz 和cumtrapz分别计算结果都与矩形公式(9.3),(9.4)计算结果的平均数z相等,即z3=z3c=z=0.395 4,梯形公式计算的结果与精确值Fs = 0.396 1更接近.但是,函数cumtrapz计算的结果展示了每次累加计算过程的结果.(四) 自编梯形数值积分的MATLAB主程序另外,也可以根据梯形公式和连续增加的
20、子区间数来逼近d 第次递归在个等距点处对的值的和的方法,现提供计算数值积分的程序如下:梯形数值积分的MATLAB主程序输入量:fun是被积函数, a,b分别是积分下、上限,是递归的次数.输出量:T是利用梯形公式数值计算d的递归值.function T=rctrap(fun,a,b,m)n=1;h=b-a; T=zeros(1,m+1); x=a; T(1)=h*(feval(fun,a)+feval(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(fun,x);endT(i+1)=T(i)/2
21、+h*s;endT=T(1:m);例9.3.6 用rctrap计算ed,递归14次,并将计算结果与精确值比较.解 (1)建立fun.m文件命名的M文件函数function y = fun(x) y = exp( (-x.2)./2)./(sqrt(2*pi);(2)输入程序>>T=rctrap(fun,0,pi/2,14), syms t fi=int(exp(-t2)/2)/(sqrt(2*pi),t,0, pi/2); Fs= double(fi), wT= double(abs(fi-T)运行后屏幕显示精确值Fs,用rctrap计算的递归值T和T与精确值Fs的绝对误差wT如下
22、T = Columns 1 through 4 0.40457385587044 0.43245899179539 0.43953669400491 0.44129851884975 Columns 5 through 8 0.44173842995173 0.44184837274713 0.44187585624611 0.44188272698315 Columns 9 through 12 0.44188444465881 0.44188487407718 0.44188498143174 0.44188500827038 Columns 13 through 14 0.4418850
23、1498004 0.44188501665745Fs = 0.44188501721659wT = Columns 1 through 4 0.03731116134615 0.00942602542121 0.00234832321168 0.00058649836684 Columns 5 through 8 0.00014658726486 0.00003664446946 0.00000916097048 0.00000229023344 Columns 9 through 12 0.00000057255779 0.00000014313941 0.00000003578485 0.
24、00000000894621 Columns 13 through 14 0.00000000223655 0.00000000055914由此可见,递归14次,用rctrap计算数值积分的结果与精确值的绝对误差为,0.441 885 017.9.3.3 辛普森(Simpson)公式及其误差分析定理9.2 (复合辛普森公式)设函数在闭区间上的个等距节点处有定义,且其函数值为 ,则在上有, (9.11)使得 d, (9.12)其中是的截断误差,即余项.(9.11)式称为复合辛普森(Simpson)数值积分公式,简称辛普森(Simpson)公式.例9.3.7 用辛普森(Simpson)公式(9.1
25、1)计算ed,取个等距节点,并将计算结果与精确值比较,然后再取计算,观察对误差的影响.解 由,得.根据辛普森(Simpson)公式(9.11)编写并输入下面的程序>> a=0;b=1;m=10000; h=(b-a)/(2*m); x=a:h:b; y=exp(-x.2)./2)./(sqrt(2*pi);z1=y(1)+y(2*m+1); z2=2*sum(y(2:2:2*m); z3=4*sum(y(3:2:2*m);z=(z1+z2+z3)*h/3, syms t,f=exp(-t2)/2)/(sqrt(2*pi);intf=int(f,t,a,b), Fs=double(i
26、ntf); Juewucha=abs(z-Fs)运行后屏幕显示用辛普森公式(9.11)计算定积分的近似值z和精确值intf及其绝对误差Juewucha(取个等距节点)如下z = 0.34133406408431intf =1125899906842624/5644425081792261*erf(1/2*2(1/2)*2(1/2)*pi(1/2)Juewucha = 1.068198423692657e-005然后再取计算的近似值z和与精确值intf及其绝对误差Juewucha的结果为:z = 0.323 261 353 494 30, Fs = 0.341 344 746 068 54, J
27、uewucha = 0.018 083 392 574 25.由此可见,节点的个数越大,步长越小,其绝对误差也越小.这个结论我们可以用证明推论9.1的方法证明,其结果见下面的推论.例9.3.8 估计用辛普森公式计算定积分ed时的误差,取.解 根据估计误差公式(9.14),先输入求的程序>>syms x,y=exp(sin(x); yx4=diff(y,x,4)运行后输出被积函数的四阶导函数. 然后在输入误差估计程序>>h=pi/40; x=0:0.00001:pi/2; yx4=sin(x).*exp(sin(x)-4*cos(x).2.*exp(sin(x)+3*si
28、n(x).2.*exp(sin(x)-6*sin(x).*cos(x).2.*exp(sin(x)+cos(x).4.*exp(sin(x);juyx4= abs(yx4); RS=(h4)*(pi/2)*max(juyx4)/180运行后屏幕显示误差估计值RS = 3.610450295892220e-0069.3.4 辛普森(Simpson)数值积分的MATLAB程序用MATLAB软件和辛普森(Simpson)公式计算数值积分和估计误差限除了用例9.3.7和例9.3.8自己编写的MATLAB计算程序外,MATLAB系统还提供了用自适应辛普森公式计算数值积分的程序quad.m,具体的调用格式
29、如下:调用格式一:quad(fun,a,b) 计算函数Y = fun(X)在区间(a,b)上的数值积分,自动选择步长,在MATLAB 6.3中绝对误差限为,在MATLAB 5.3中误差限为,输出数值积分值.其中函数Y = fun(X)是以fun.m文件命名的M文件函数(或库函数如sin,log),或者是F = inline(fun)形式,或者数组,函数Y = fun(X)将接受向量的自变量,并且返回结果是向量Y,即在X的每个元素处的积分估计值.调用格式二:quad(fun,a,b,tol) 同上,但指定绝对误差限为tol,默认值为 ,在MATLAB 5.3中误差限为.调用格式三:Q,FCNT
30、= quad (.) 输出数值积分值Q和函数赋值的编号FCNT.调用格式四:quad(fun,a,b, tol,TRACE)输入量非零数TRACE,则会以动态的形式显示在递归求积分的整个过程的fcnt a b-a Q的值; 调用格式五:quad(fun,a,b, tol,TRACE,P1,P2, )此命令规定对函数fun (X,P1,P2,) 附加的参数P1, P2, .直接通过.传递tol 或 TRACE的空矩阵使用默认值.我们还可以根据复合辛普森(Simpson)公式编写名为comsimpson.m的计算d 的数值积分的MATLAB主程序.复合辛普森(Simpson)数值积分的MATLAB
31、主程序输入量:fun是被积函数, a,b分别是积分下、上限,是小区间的数目.输出量:y是利用复合辛普森数值积分公式(9.11)计算d的数值积分结果.function y=comsimpson(fun,a,b,n) z1=feval (fun,a)+ feval (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 (fun,x2); endfor k=3:2:2*m x3=x+k*h; z3= z3+4*feval (fun,x3); endy=(z1+z2+
32、z3)*h/3;例9.3.9 用comsimpson.m和quad.m分别计算定积分ed,取精度为,并与精确值比较.解 输入程序>> Q1,FCNT14 = quad(fun,0,1,1.e-4,3), Q2 =comsimpson (fun,0,1,10000)syms x fi=int(exp( (-x.2)./2)./(sqrt(2*pi),x,0, 1); Fs= double (fi)wQ1= double (abs(fi-Q1) ), wQ2= double (abs(fi-Q2) )运行后屏幕显示的精确值Fs,用comsimpson.m和quad.m分别计算的近似值Q
33、2、Q1和迭代次数FCNT14,取精度分别为,Q2、Q1分别与精确值Fs的绝对误差wQ2, wQ1如下9 0.0000000000 2.71580000e-001 0.1070275100 11 0.2715800000 4.56840000e-001 0.1597942242 13 0.7284200000 2.71580000e-001 0.0745230082Q1 = FCNT14 = Q2 = 0.3413 13 0.3413Fs = wQ1 = wQ2 = 0.3413 3.6619e-009 3.7061e-0059.3.6 牛顿-科茨(Newton-Cotes)数值积分和误差分析
34、的MATLAB程序下面介绍计算科茨(Cotes)系数和估计(9.21)式在上的截断误差的程序及其用法,并且介绍根据在MATLAB中提供了用自适应递归牛顿科茨公式计算数值积分的程序quad8.m和使用方法.(一) 估计误差的MATLAB程序根据计算定积分d的近似值的n阶牛顿-科茨公式的截断误差公式(9.24)式和(9.25)式,现提供的求(9.21)式在上的截断误差公式的主程序如下:计算n阶牛顿-科茨的公式的截断误差公式的MATLAB主程序输入量: 牛顿-科茨的公式(9.21)的阶数.输出量: RNC是对应的n阶牛顿-科茨的公式(9.21)的截断误差公式.function RNC=ncE(n)s
35、uk=1; p=n/2-fix(n/2);if p=0for k=1:n+2suk=suk*k;endsuk; syms t a b fxn2,su=t2; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC= (b-a)/n)(n+3)*fxn2*abs(y)/ suk;elsefor k=1:n+1suk=suk*k;endsuk; syms t a b fxn1,su=t; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC=
36、 (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*fxn1n = RNC2 = 2 1/90*(1/2*b-1/2*a)5*fxn2n = RNC3 = 3 3/8
37、0*(1/3*b-1/3*a)5*fxn1n = RNC4 = 4 8/945*(1/4*b-1/4*a)7*fxn2n = 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 xy=exp(sin(x); yx6=diff(y,x,6), yx8=dif
38、f(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)-
39、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
40、(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*myx6RNC6 =9/1400*(1/6*b-1/6*a)9*m
41、yx8运行后屏幕显示误差限如下 RNC5 = RNC6 = 3.625385713339320e-004 3.826183594914823e-005(二) 计算科茨(Cotes)系数和求截断误差公式的MATLAB程序现提供计算定积分I=d的n阶牛顿科茨的公式(9.21)、(9.23)编写的计算科茨(Cotes)系数主程序如下:计算n阶科茨(Cotes)系数和求截断误差公式的MATLAB主程序输入量:牛顿科茨的公式(9.21)的阶数.输出量:Cn是计算定积分I的n阶牛顿科茨的公式的科茨(Cotes)系数,RNCn是对应的n阶牛顿科茨的公式(9.21)的截断误差公式.function Cn, R
42、NCn=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:nsu1=su*(t-u);m01=m0*u; su=su1;m0=m01;endsu;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:nk1=k*j; m1=m*(n-j); f=su/(t-j); intf=int(f,t,0,n); y=double(intf);Cn(
43、j+1)=(-1)(n-j)*y)/(n*k1*m1);warning off MATLAB:divideByZeroend 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=0for k=1:n+2suk=suk*k;endsuk; syms t a b fxn2,su=t2; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNCn=(b-a)/n)(n+3)*fxn
44、2*abs(y)/suk;elsefor k=1:n+1suk=suk*k;endsuk; syms t a b fxn1,su=t; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNCn=(b-a)/n)(n+2)*fxn1*abs(y)/ suk;end例9.3.16 用计算n阶科茨(Cotes)系数和求截断误差公式的MATLAB主程序,计算定积分I=d的阶牛顿科茨公式的系数和截断误差公式.解 (1)先求阶牛顿科茨公式的系数和截断误差公式.输入程序>> n1=1,Cn1,RNCn1=newcotE
45、(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*fxn1n2 = 2 Cn2 = 1/6 2/3 1/6 RNCn2 =1/90*(1/2*b-1/2*a)5*fxn2n3 = 3 Cn3 = 1/8 3/8 3/8 1/8 RNCn3 =3/80*(1/3*b-1/3*a)5*fxn1(三) 计算牛顿科茨公式的MATLAB程序在MATL
46、AB5.3中提供了用自适应递归牛顿科茨公式计算数值积分的程序quad8.m.但是,因为quad8.m递归速度慢和其它原因,在MATLAB6.5和MATLAB7.01中极力推荐使用程序quadl.m取代quad8.m. 在MATLAB7.01中虽然还保留quad8.m的文件名,实际上已经将它的程序废弃,输入有关quad8的程序时,运行的是quadl.m的程序.quad8.m和quadl.m的调用方法与quad 类似,下面介绍quad8.m的调用格式.调用格式一:quad8(fun,a,b) 计算函数Y = fun(X)在区间 (a,b)上的数值积分,自动选择步长,相对误差限为,输出数值积分值.函
47、数Y = fun(X)将接受向量的自变量X,并且返回结果是向量Y,即在X的每个元素处的积分估计值.调用格式二:quad8(fun,a,b,tol) 同上,返回相对误差tol的数值积分,如果TOL = rel_tol abs_tol,则rel_tol 和abs_tol 分别表示相对误差和绝对误差.调用格式三: quad8(fun,a,b, tol,TRACE)输入非零数TRACE,则会以动态的形式显示误差为tol的数值积分在递归求积分的整个过程和图形. 调用格式四:quad8(fun,a,b, tol,TRACE,P1,P2,.)此命令规定对函数fun (X,P1,P2,.) 附加的参数P1,
48、P2, .直接通过.传递tol 或 TRACE的空矩阵使用默认值.例9.3.17 用梯形公式、辛普森(Simpson)和牛顿科茨求积公式计算定积分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.10437572798742ztc = Columns
49、 1 through 3 0 0.00630298652792 0.01264569951380.Columns 250 through 251 3.08729642810745 3.10437572798742(2)用辛普森(Simpson)求积公式计算定积分. 输入程序>> syms xL= inline(' exp(sin(x)'); QS,FCNTS =quad(L,0, pi/2,1.e-4,2)运行后屏幕显示用辛普森(Simpson)求积公式计算定积分的值QS和递归次数FCNTS分别如下QS = FCNTS = 3.10438133817254 13(3
50、)用牛顿科茨求积公式计算定积分. 在MATLAB6.5中输入程序>> syms xL= inline(' exp(sin(x)'); Q8,FCNT8 = quad8(L,0, pi/2,1.e-4,3)运行后屏幕显示用牛顿科茨求积公式计算定积分的值Q8和递归次数FCNTS分别如下 Q8 = FCNT8 = 3.10437901785555 33(4)输入求定积分的精确值的程序>> syms xy=exp(sin(x); F=int(y,0,pi/2), Fj=double(F)wzt=abs( Fj- zt), wQS = abs(Fj- QS), w
51、Q8 = abs(Fj- Q8)运行后屏幕显示计算的定积分值F和近似值Fj,梯形公式、辛普森(Simpson)和牛顿科茨求积公式计算定积分的值与Fj的绝对误差wztc, wQS和wQ8如下Warning: Explicit integral could not be found.> In C:MATLAB6p5p1toolboxsymbolicsymint.m at line 58F =int(exp(sin(x),x = 0 . 1/2*pi)Fj = wzt = 3.10437901785556 3.289868133471430e-006wQS = wQ8 = 2.32031698
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- (3篇)五年级家长会发言稿格式范文
- 运动会跳高广播稿15篇
- 销售内勤工作总结15篇
- 连锁业务工作总结
- 送清凉活动工作总结
- 中华民族发展史知到课后答案智慧树章节测试答案2025年春云南大学
- 小学作文课件:如何写好关于西瓜的作文
- 人教山西 九年级 下册 语文 第二单元《 山西中考题型专练》习题课 课件
- 人教版高中语文第二册《物种起源》导言 同步练习基础部分
- 高中语文必修3巩乃斯的马 同步练习3
- 2025年广东省佛山市南海区中考一模英语试题(原卷版+解析版)
- 2025河北石家庄市国有企业招聘21人笔试参考题库附带答案详解
- 单独招生机电类试题库含答案
- 国开2025年春季《形势与政策》大作业答案
- 上海市农村房地一体宅基地确权登记工作实施方案
- 计算机网络知到智慧树章节测试课后答案2024年秋贵州财经大学
- 《行业分析方法》课件
- 屋面光伏工程施工组织设计
- 小学校园欺凌治理委员会
- 互联网护理服务典型案例
- Unit 3 Keep fit 知识点课件 合作探究一
评论
0/150
提交评论