Matlab中上下标及希腊字母.docVIP

很多时候都要在matlab画图的时候添加一些公式符号之类的，有一些特殊的字符并不能直接从键盘上输入，比如希腊字母等等。但是有想用，因为这样使图看起来漂亮而且容易理解。 例如：我想输入摄氏度的符号，怎么办咧？ 也许你突然想到，摄氏度不就是一个小圆圈加一个大写的C么。 于是就用T=25oC来表示了，一看，多少还算是那么回事，但怎么看怎么有点别扭。 因为o作为上标的时候它不是一个正真的圆圈，最多是个椭圆，并且它体积太大了。 好吧，既然如此，那就用中文输入法打个句号“。”上去，即T=25。C 再看的时候，圆圈到是圆圈了，但还是别扭呢？ 那是因为在编码中，中文句号占了两个字符的位置，所以圆圈和C的位置拉得太远，根本看不下去。 行了，告诉你吧 正确的表示方法为：T=25circC，这样就好看多了！下面给出Matlab中下标及希腊字母的使用方法，还有更多的使用方法可以参考matlab帮助 文档中的Text Properties: 下标用 _(下划线) 上标用 （尖号） 希腊字母等特殊字符用 加拼音如如果是需要大写的希腊字符，加相应的名称首字母改为大写。比如： theta Theta kappa rho alpha beta gamma theta Theta Gamma delta Delta xi Xi elta epsilong zeta miu nu tau lamda Lamda pi Pi sigma Sigma phi Phi psi Psi chi omega Omega geq 不等于 neq gg 正负 pm 左箭头 leftarrow 右箭头 rightarrow 上箭头 uparrow 上圆圈（度数） circ 例 text(2,3,alpha_2beta) 注： 可用把须放在一起的括起来MATLAB基本绘图函数plot: x轴和y轴均为线性刻度（Linear scale）loglog: x轴和y轴均为对数刻度（Logarithmic scale）semilogx: x轴为对数刻度，y轴为线性刻度semilogy: x轴为线性刻度，y轴为对数刻度=若要画出多条曲线，只需将座标对依次放入plot函数即可：plot(x, sin(x), x, cos(x);若要改变颜色，在座标对後面加上相关字串即可：plot(x, sin(x), c, x, cos(x), g);若要同时改变颜色及图线型态（Line style），也是在座标对後面加上相关字串即可：plot(x, sin(x), co, x, cos(x), g*);=小整理：plot绘图函数的叁数字元 颜色 字元 图线型态y 黄色 . 点k 黑色 o 圆w 白色 x xb 蓝色 + +g 绿色 * *r 红色 - 实线c 亮青色 : 点线m 锰紫色 -. 点虚线 - 虚线=图形完成後，我们可用axis(xmin,xmax,ymin,ymax)函数来调整图轴的范围：axis(0, 6, -1.2, 1.2);此外，MATLAB也可对图形加上各种注解与处理：xlabel(Input Value); % x轴注解ylabel(Function Value); % y轴注解title(Two Trigonometric Functions); % 图形标题legend(y = sin(x),y = cos(x); % 图形注解grid on; % 显示格线我们可用subplot来同时画出数个小图形於同一个视窗之中：subplot(2,2,1); plot(x, sin(x);subplot(2,2,2); plot(x, cos(x);subplot(2,2,3); plot(x, sinh(x);subplot(2,2,4); plot(x, cosh(x);MATLAB还有其他各种二维绘图函数，以适合不同的应用，详见下表。=小整理：其他各种二维绘图函数bar 长条图errorbar 图形加上误差范围fplot 较精确的函数图形polar 极座标图hist 累计图rose 极座标累计图stairs 阶梯图stem 针状图fill 实心图feather 羽毛图compass 罗盘图quiver 向量场图=以下我们针对每个函数举例。当资料点数量不多时，长条图是很适合的表示方式：close all; % 关闭所有的图形视窗x=1:10;y=rand(size(x);bar(x,y);如果已知资料的误差量，就可用errorbar来表示。下例以单位标准差来做资料的误差量：x = linspace(0,2*pi,30);y = sin(x);e = std(y)*ones(size(x);errorbar(x,y,e)对於变化剧烈的函数，可用fplot来进行较精确的绘图，会对剧烈变化处进行较密集的取样，如下例：fplot(sin(1/x), 0.02 0.2); % 0.02 0.2是绘图范围若要产生极座标图形，可用polar：theta=linspace(0, 2*pi);r=cos(4*theta);polar(theta, r);对於大量的资料，我们可用hist来显示资料的分情况和统计特性。下面几个命令可用来验证randn产生的高斯乱数分：x=randn(5000, 1); % 产生5000个 ?=0，?=1 的高斯乱数hist(x,20); % 20代表长条的个数rose和hist很接近，只不过是将资料大小视为角度，资料个数视为距离，?眉昊嬷票硎荆?x=randn(1000, 1);rose(x);stairs可画出阶梯图：x=linspace(0,10,50);y=sin(x).*exp(-x/3);stairs(x,y);stems可产生针状图，常被用来绘制数位讯号：x=linspace(0,10,50);y=sin(x).*exp(-x/3);stem(x,y);stairs将资料点视为多边行顶点，并将此多边行涂上颜色：x=linspace(0,10,50);y=sin(x).*exp(-x/3);fill(x,y,b); % b为蓝色feather将每一个资料点视复数，并以箭号画出：theta=linspace(0, 2*pi, 20);z = cos(theta)+i*sin(theta);feather(z);compass和feather很接近，只是每个箭号的起点都在圆点：theta=linspace(0, 2*pi, 20);z = cos(theta)+i*sin(theta);compass(z);-3.基本XYZ立体绘图命令在科学目视表示（Scientific visualization）中，三度空间的立体图是一个非常重要的技巧。本章将介绍MATLAB基本XYZ三度空间的各项绘图命令。mesh和plot是三度空间立体绘图的基本命令，mesh可画出立体网状图，plot则可画出立体曲面图，两者产生的图形都会依高度而有不同颜色。下列命令可画出由函数 形成的立体网状图:x=linspace(-2, 2, 25); % 在x轴上取25点y=linspace(-2, 2, 25); % 在y轴上取25点xx,yy=meshgrid(x, y); % xx和yy都是21x21的矩阵zz=xx.*exp(-xx.2-yy.2); % 计算函数值，zz也是21x21的矩阵mesh(xx, yy, zz); % 画出立体网状图surf和mesh的用法类似：x=linspace(-2, 2, 25); % 在x轴上取25点y=linspace(-2, 2, 25); % 在y轴上取25点xx,yy=meshgrid(x, y); % xx和yy都是21x21的矩阵zz=xx.*exp(-xx.2-yy.2); % 计算函数值，zz也是21x21的矩阵surf(xx, yy, zz); % 画出立体曲面图为了方便测试立体绘图，MATLAB提供了一个peaks函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点，其方程式为：要画出此函数的最快方法即是直接键入peaks：peaksz = 3*(1-x).2.*exp(-(x.2) - (y+1).2) .- 10*(x/5 - x.3 - y.5).*exp(-x.2-y.2) .- 1/3*exp(-(x+1).2 - y.2)我们亦可对peaks函数取点，再以各种不同方法进行绘图。meshz可将曲面加上围裙：x,y,z=peaks;meshz(x,y,z);axis(-inf inf -inf inf -inf inf);waterfall可在x方向或y方向产生水流效果：x,y,z=peaks;waterfall(x,y,z);axis(-inf inf -inf inf -inf inf);下列命令产生在y方向的水流效果：x,y,z=peaks;waterfall(x,y,z);axis(-inf inf -inf inf -inf inf);meshc同时画出网状图与等高线：x,y,z=peaks;meshc(x,y,z);axis(-inf inf -inf inf -inf inf);surfc同时画出曲面图与等高线：x,y,z=peaks;surfc(x,y,z);axis(-inf inf -inf inf -inf inf);contour3画出曲面在三度空间中的等高线：contour3(peaks, 20);axis(-inf inf -inf inf -inf inf);contour画出曲面等高线在XY平面的投影：contour(peaks, 20);plot3可画出三度空间中的曲线：t=linspace(0,20*pi, 501);plot3(t.*sin(t), t.*cos(t), t);亦可同时画出两条三度空间中的曲线：t=linspace(0, 10*pi, 501);plot3(t.*sin(t), t.*cos(t), t, t.*sin(t), t.*cos(t), -t);y(2:4)-1 % 取出y的第二至第四个元素来做运算同样的方法可用於产生公差为1的等差数列：x = 7:16x = 7:3:16 % 公差为3的等差数列x = linspace(4, 10, 6) % 等差数列：首项为4,末项为10,项数为6若要重新安排矩阵的形状，可用reshape命令：B = reshape(A, 4, 2) % 4是新矩阵的列数，2是新矩阵的行数举例来说，下列命令会产生一个长度为6的调和数列（Harmonicsequence）：x = zeros(1,6); % x是一个16的零矩阵for i = 1:6,x(i) = 1/i;endfor圈可以是多层的，下例产生一个16的Hilbert矩阵h，其中为於第i列、第j行的元素为：h = zeros(6);for i = 1:6for j = 1:6h(i,j) = 1/(i+j-1);endendformat rat % 使用分数来表示数值disp(x)1 1/2 1/3 1/4 1/5 1/6function output = fact(n)% FACT Calculate factorial of a given positive integer.output = 1;for i = 1:n,output = output*i;end其中fact是函数名，n是输入引数，output是输出引数，而i则是此函数用到的暂时变数。要使用此函数，直接键入函数名及适当输入引数值即可：MATLAB的函数也可以是递式的（Recursive），也就是说，一个函数可以呼叫它本身。举例来说，n! =n*(n-1)!，因此前面的阶乘函数可以改成递式的写法：function output = fact(n)% FACT Calculate factorial of a given positive integer recursively.if n = 1, % Terminating conditionoutput = 1;return;endoutput = n*fact(n-1);在写一个递函数时，一定要包含结束条件（Terminatingcondition），否则此函数将会一再呼叫自己，永远不会停止，直到电脑的记忆体被耗尽为止。以上例而言，n=1即满足结束条件，此时我们直接将output设为1，而不再呼叫此函数本身。勒让德函数/view/ab40bb4cf7ec4afe04a1df81.html第九章 数值积分9.1 积分的MATLAB符号计算9.1.1 定积分的MATLAB符号计算9.1.5 由，所围成的平面区域D.求平面区域D的面积S.解 输入作函数图形的程序 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), ylabel(y),title(y=sinx , y=cosx 和x=-0.5及x=1.5所围成的平面区域的图形)运行后屏幕显示图形.求平面区域D的面积S.输入计算面积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 及其近似值 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，求.解 输入程序： 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)运行后屏幕显示计算变限积分的导数.9.2 数值积分的思想及其MATLAB程序9.2.3 矩形公式的MATLAB程序（一） 函数sum的调用格式调用格式一：sum(X)调用格式二：sum (X,DIM)例9.2.2 用MATLAB和矩形公式（9.3）、（9.4）计算ed，并与精确值比较.解 将分成20等份，步长为，输入程序 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求定积分的精确值，输入程序 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）计算结果的绝对误差wz1、wz2.（二） 函数cumsum的调用格式调用格式一：cumsum(X) 调用格式二：cumsum (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(20)*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.3961wz1 = wz2 = 0.0088 0.00759.3 插值型数值积分及其MATLAB 程序9.3.2 梯形公式的MATLAB程序(一) 根据梯形公式和估计误差公式自己编写MATLAB程序计算定积分例9.3.2 分别取，用梯形公式计算定积分ed，并与精确值比较.然后观察对计算结果的有效数字和绝对误差的影响.解 编写并输入如下程序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(z8000-Fs)运行后屏幕显示取时，积分区间上等距节点的个数，用梯形公式计算定积分的值z8000和精确值intf的近似值Fs及其绝对误差Juewucha8000. (二) 用函数trapz计算定积分调用格式一：Z =trapz(Y)调用格式二：Z =trapz(X,Y)调用格式三：Z = trapz (X,Y,DIM) 或 trapz (Y,DIM) (三) 用函数cumtrapz计算定积分调用格式一：Z =cumtrapz (Y)调用格式二：Z =cumtrapz (X,Y)调用格式三：Z = cumtrapz (X,Y,DIM) 或 cumtrapz (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)*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.(四)梯形数值积分的MATLAB主程序梯形数值积分的MATLAB主程序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+h*s;endT=T(1:m);例9.3.6 用rctrap计算ed，递归14次，并将计算结果与精确值比较.解 输入程序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的绝对误差wT9.3.3 辛普森公式及其误差分析例9.3.7 用辛普森公式计算ed，取个等距节点，并将计算结果与精确值比较，然后再取计算，观察对误差的影响.解 由，得.根据辛普森（Simpson）公式编写并输入下面的程序 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(intf); Juewucha=abs(z-Fs)运行后屏幕显示用辛普森公式(9.11)计算定积分的近似值z和精确值intf及其绝对误差Juewucha（取个等距节点）.例9.3.8 估计用辛普森公式计算定积分ed时的误差，取.解 根据估计误差公式，先输入求的程序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*sin(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程序调用格式一：quad(fun,a,b)调用格式二：quad(fun,a,b,tol)调用格式三：Q,FCNT = quad (.) 调用格式四：quad(fun,a,b, tol,TRACE)调用格式五：quad(fun,a,b, tol,TRACE,P1,P2, )复合辛普森（Simpson）数值积分的MATLAB主程序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+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分别计算的近似值Q2、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）数值积分和误差分析的MATLAB程序(一) 估计误差的MATLAB程序计算n阶牛顿-科茨的公式的截断误差公式的MATLAB主程序function RNC=ncE(n)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);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= (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/80*(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=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*myx6RNC6 =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: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(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)*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);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*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程序调用格式一：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.10437572798742ztc = Columns 1 through 3 0 0.00630298652792 0.01264569951380.Columns 250 through 251 3.08729642810745 3.10437572798742（2）用辛普森求积公式计算定积分. 输入程序 syms xL= 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 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), wQ8 = abs(Fj- Q8)运行后屏幕显示计算的定积分值F和近似值Fj，梯形公式、辛普森和牛顿科茨求积公式计算定积分的值与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.320316987436399e-006 7.993605777301127e-015（5）输入画图程序 syms xF=int(exp(sin(x),0,pi/2),Fj=double(F),plot(Fj,ro), hold onL= 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 offhold on, Q8,FCNT8 = quad8(L,0, pi/2,1.e-4,3), hold offgrid, xlabel(自变量 X), ylabel(因变量 Y)title(牛顿科茨公式和梯形公式计算