微分方程的数值积分课件.ppt

1、5.1 函数极限和导数函数极限和导数5.1.1 单变量函数值的计算和绘图单变量函数值的计算和绘图【例例5-1-1】设要求以0.01s为间隔，求出y的151个点，并求出其导数的值和曲线。解解:建模可以采取下列两种方法来做:(1)只用主程序编程的方法;(2)编成函数文件，由主程序调用的方法。求导数可采用diff函数对数组y做运算的方法。3t 34sine23y4t MATLAB程序(1)第一种方法的程序exn511a如下:t=0:.01:1.5;%设定自变量数组t w=4*sqrt(3);%固定频率 y=w/8*exp(-4*t).*sin(w*t+pi/3);%注意用数组运算式 subplot(

2、2，1，1)，plot(t，y)，grid%绘制曲线并加上坐标网格 title(绘图示例)，xlabel(时间 t)，ylabel(y(t)%加标注 Dy=diff(y);subplot(2，1，2)，plot(t(length(t)-1)，Dy)，grid%求导数并绘制曲线，注意用数值方法求导后，导数数组长度比原函数减少一 ylabel(Dy(t)%加标注(2)第二种方法的主程序exn511b如下:dt=0.01;t=0:dt:1.5;w=4*sqrt(3);y=exn511bf(t，w);Dy=diff(y)/dt;%绘图和加标注的程序略去 另要建立一个函数文件exn511bf.m，其内容

3、为:function xvalues=exn511bf(tvalues，w)%注意编写的函数文件中，其语法对数组w应该能用元素群运算。xvalues=w*exp(-4*tvalues).*sin(w*tvalues+pi/3);程序运行结果运行这两种程序都得到图5-1所示的曲线。为了节省篇幅，我们没有显示y的数据。以后的各例中还将省略绘图时的标注语句。从本例看，第二种方法似乎更麻烦一些，但它具备模块化的特点。当程序中要反复多次调用此函数，而且输入不同的自变量时，利用函数文件可大大简化编程。我们应该掌握这种方法。两次应用diff函数或用diff(y，2)可以求y的二次导数，读者可自行实践。图 5

4、-1 例5-1-1的曲线5.1.2 参变方程表示的函数的计算和绘图参变方程表示的函数的计算和绘图【例例5-1-2】摆线的绘制。如图5-2所示，当圆轮在平面上滚动时，轮上任一点所画出的轨迹称为摆线。如果这一点不在圆周上而在圆内，则生成内摆线;如果该点在圆外，即离圆心距离大于半径，则生成外摆线。对于后一种情况，可由火车轮来想象。其接触轨道的部分，并不是直径最大处，其内侧的直径还要大一些，以防止车轮左右出轨。在这部分边缘上的点就形成外摆线。图 5-2 摆线的生成解：解：建模概括几种情况，其普遍方程可表示为:xA=rt-Rsin t yA=r-Rcos t可由这组以t为参数的方程分析其轨迹。MATLA

5、B程序 t=0:0.1:10;%设定参数数组 r=input(r=)，R=input(R=)%输入常数 x=r*t-R*sin(t);%计算x，y y=r-R*cos(t);plot(x，y)，axis(equal)%绘图 程序运行结果设r=1，令R=r，R=0.7及R=1.5时得到的摆线、内摆线和外摆线都绘于图5-3中。为了显示摆线的正确形状，x、y坐标保持等比例是很重要的，因此程序中要加axis(equal)语句。图 5-3 摆线的绘制5.1.3 曲线族的绘制曲线族的绘制【例例5-1-3】三次曲线的方程为y=ax3+cx，试探讨参数a和c对其图形的影响。解解:建模因为函数比较简单，因此可以

6、直接写入绘图语句中，用循环语句来改变参数。注意坐标的设定方法，以得到适于观察的图形。给出的程序不是唯一的，例如也可用fplot函数等。读者可自行探索其他编法。MATLAB程序 x=-2:0.1:2;%给定x数组，确定范围及取点密度 subplot(1，2，1)%分两个画面绘图 for c=-3:3%取不同的c循环 plot(x，x.3+c*x)，hold on，end，grid subplot(1，2，2)for a=-3:3%取不同的a循环 plot(x，a*x.3+x)，hold on，end，grid，hold off 用直接在图形窗内编辑的方法可在图内标注字符。程序运行结果程序运行结果

7、见图5-4，其中a和c均从-3取到3，步长为1。图 5-4 c和a取不同值时y=ax3+cx的曲线族(a)a=1，c取不同值;(b)c=1，a取不同值5.1.4 极限判别极限判别用MATLAB来表达推理过程是比较困难的，因为它必须与实际的数值联系起来。比如无法用无穷小和高阶无穷小的概念，只能用e-10、e-20等数值。不过极限的定义恰恰用了和这些数的概念，因此不难用程序表达。【例例5-1-4】极限的定义和判别。解解:建模 对于函数y=f(x)，当任意给定一个正数时，有一个对应的正数存在，使得当0|xc-x|epsilon delta=delta/2;x=xc-delta;%找 if abs(d

8、elta)0.0001 f=x0.3+10*x0.2-2*sin(x0)-50;%求f(x0)g=3*x02+20*x0-2*cos(x0);%求函数的导数g(x0)xx=x0-f/g;%切线法求解公式 e=abs(xx-x0);x0=xx%精度控制,把新值赋予x0 plot(x0，0，*)，pause(2)%画出新点在x轴上的位置 end 程序运行结果运行此程序，设x0分别为-10、-2、2，得出的解分别为-9.4384、-2.5648、2.0707,从图上可看到解逐渐接近精确解的过程。本程序中用了导数的解析形式，对于某些较为复杂的函数,这样会很繁琐，所以这不是一个值得推荐的程序，只是为了帮

9、助读者理解，好的程序应该用数值方法求导数(如例5-1-1)，或用两分法来求根。读者可自行改进这个程序。懂得原理后，实际上不必再去编程，直接用fzero求f(x)的根即可，要做的工作是把f(x)定义为一个用M文件表示的函数，例如建立一个exn517f.m文件，其内容为function y=exn517f(x)y=x.3+10*x.2-2*sin(x)-50;再调用x=fzero(exn517f，x0)，即可得出解x。【例例5-1-8】设f(x)=x3+2x，x0=1。(a)在x0-1/2xx0+3的区间内画出y=f(x)的曲线;(b)保持x0不动，求差商 ，作为h的函数;(c)求h0时q(h)的

10、极限;(d)对于h=3、2、1，定义割线y=f(x0)+q*(x-x0)，画出此割线;(e)在x0=0、1、2三个点上分别画出曲线的切线。h)f(xh)f(xq(h)00解:建模先在指定的定义域内画出函数的图形，然后在x0附近计算不同h时的差商及其在h0时的极限，再按切线公式画出切线图形。MATLAB程序 x=linspace(-1/2，3，100);subplot(1，2，1)，plot(x，x.3+2*x)，grid on，hold on%画曲线在一个点上，取不同步长求斜率，画出多根割线 h=1，0.1，0.01;for k=1:3 x1=1，1+h(k);f=x1.3+2*x1;%取不同

11、步长 q=diff(f)/diff(x1)，%求差商 plot(x，f(1)+q*(x-x1(1)，:)，%画割线 end%在多个点上，取同样步长h=0.1，求斜率，画出多根切线 for k=0:2 x1=k，k+0.1;f=x1.3+2*x1;%取不同位置的x，求函数 q=diff(f)/diff(x1)，%求差商 subplot(1，2，2)，plot(x，x.3+2*x)，grid on，hold on plot(x，f(1)+q*(x-x1(1)，:r)，%画割线 end 程序运行结果在x0=1处取三个不同的步长h，得到的斜率分别为:q=9 5.3100 5.0301用这三个斜率画出的

12、近似切线如图5-10(a)所示，可见q逐步接近精确值。取三个不同的x0，得到的斜率分别为:q=2.0100 5.3100 14.6100根据这些斜率画出的切线如图5-10(b)所示。图 5-10 曲线的割线与切线的关系斜率的精确值只能通过符号数学来求，其程序为syms x，y=x3+2*x，Dy=diff(y)，Dy1=subs(Dy，x，0，1，2)运行此程序的结果为Dy1=2 5 14【例例5-1-9】求函数 位于区间0，内的近似极值。解解:先完全采用数值方法来解这个问题，为了尽量求得准确一些，对自变量作较细的分度，求出其函数值，画图并求其极值。编写程序如下:clear，clf，x=lin

13、space(0，pi，10000);y=2*sin(2*x).2+5/2*x.*cos(x/2).2;plot(x，y)，grid on，hold on2xxcos25(2x)2siny22程序运行后得到图5-11所示的曲线。求y的最大值可有几种方法:(1)方法一:由y的数据直接求极值(exn519a):ymax1=max(y)求得 ymax=3.7323如果要进一步确定此处的x坐标和其他几个极值，那就要麻烦一些，要利用find函数找到变量的下标:%找最大值处y的下标 nm1=find(y=max(y)%求该下标处的x的值 xm1=x(nm1)运行的结果是nm1=2752，xm1=0.8643

14、。图5-11 曲线的极值位置其他两个极值点是局部极值，只能在划分出的局部区域中寻找，所以先任取两个边界点x01=1和x02=2，该两点处的下标值也要用find函数来求:n01=find(abs(x-1)2e-4);%找两个边界处的y的下标 n02=find(abs(x-2)2e-4);ymin=min(y(n01:n02)%求x=12之间y的最小值 ymax2=max(y(n02:end)%求x=2pi之间y的最大值 nm2=find(y(n01:n02)=ymin)+n01-1%求y的最小值处的下标 nm3=find(y(n02:end)=ymax2)+n02-1%求y的最大值处的下标 xm

15、2=x(nm2)，xm3=x(nm3)%求y的最小最大值处的x值程序运行结果为：ymin=1.9446 ymax2=2.9571 nm2=5171 nm3=7147 7148 xm2=1.6244 xm3=2.2452 2.2455这种方法比较直观，也有利于思考。但键入量多，也比较麻烦。MATLAB已提供了求极小值的函数fminbnd，可以用一条语句求出极小值(参阅4.4.2节)。求极大值时只需将函数的值反号。因此键入下面三条语句就可代替前面的八条语句。xm1,mym1=fminbnd(-(2*sin(2*x).2+5/2*x.*cos(x/2).2),0,1)xm2,ym2=fminbnd(

16、2*sin(2*x).2+5/2*x.*cos(x/2).2,1,2)xm3，mym3=fminbnd(-(2*sin(2*x).2+5/2*x.*cos(x/2).2),2,3)xm=xm1，xm2，xm3，ym=-mym1，ym2，-mym3程序运行结果为:xm=0.8642 1.6239 2.2449 ym=-3.7323 1.9446 2.9571把mym1、mym3反号得到ym1、ym3，因为它们是为了利用fminbnd函数才变了号的。(2)方法二:上述方法完全没有利用极值处导数为零的性质，本方法则利用导数为零性质，但用数值方法求导，其程序如下(exn519b):h=pi/(1000

17、0-1);Dy=diff(y)/h;%先用数值方法求导数 plot(x(1:end-1)，Dy，-.r)%绘制导数曲线导数的长度比原函数的长度小1，绘图时必须相应地把x的长度减小。而求出导数过零点的下标后要加1，才能用做函数的下标检测的门限值需要经过反复试验。因为x的取值是离散的，不但不可能正好取到Dy=0的点，而且用户也不知道最靠近零的那个点的Dy值有多大，它取决于x的步长和过零点处的斜率，所以设它为某一常数kh与步长h的乘积，让kh可以人工输入。检测出的过零点太多，说明kh太大；反之，则kh太小。在本题中，理想的情况是得到三个过零点。kh=input(过零点检测门限取h的kh倍，kh=)n

18、m=find(abs(Dy)kh*h)+1，%求导数过零点 ym=y(nm)%求极值试验的结果是:若kh=1，则得出的nm只包含两项，意味着只找到两个过零点，有一个没找到，需加大;若kh=8，则nm又出现了4项，后两项实际上是一个零点。于是有:nm=find(abs(Dy)8*h)，ym=y(nm)，xm=x(nm)n=2751 5169 7145 7146 ym=3.7323 1.9446 2.9571 2.9571 xm=0.8640 1.6237 2.2446 2.2449最后两项代表同一极值，可以去掉一个，得到所需的三个极值解。(3)方法三：以下提供了用符号数学工具箱解本题的程序(ex

19、n519c)，它的好处是求过零点时不必反复试验。clear，syms x，y=2*sin(2*x)2+5/2*x*cos(x/2)2;Dy=diff(y);ezplot(y,0，pi)，hold on，ezplot(Dy,0，pi),grid on下面要根据Dy=0的条件解出相应的xm值。因为这是一个超越函数，求其零点必须用近似算法，根据图形，在各极值点附近取一起始点，逐次使用fsolve函数来求。根据图形可知y在x=1、1.5、2.2附近各有一个极值点，用三条语句分别求出:x1=fsolve(4*sin(4*x)+5/4*cos(x)+5/4-5/4*x*sin(x)，1)x2=fsolve

20、(4*sin(4*x)+5/4*cos(x)+5/4-5/4*x*sin(x)，1.5)x3=fsolve(4*sin(4*x)+5/4*cos(x)+5/4-5/4*x*sin(x)，2.2)答案分别为:0.8642 1.6239 2.2449 ym=subs(y，x，x1，x2，x3)，plot(x1，x2，x3，ym，o)，grid on得极值分别为:ym=3.7323 1.9446 2.9571【应用篇中与本节相关的例题】【例6-2-2】平面上质点运动轨迹的绘制，即知道参数方程x=x(t)及y=y(t)，画出y=f(x)。【例6-2-3】平面上质点运动轨迹及李萨如图形的绘制。与上题相仿

21、，其特点在于这两个参数方程都是周期性的函数，因此得出的是李萨如图形。【例6-6-1】声波的波形叠加造成的拍频。两个波形分别为x1=f1(t)及x2=f2(t)，两个波形叠加后的x=x1+x2=f1(t)+f2(t)，画出图形就可以直接看出其合成效果和拍频，也可使之变成可听到的声音而直接感受到。【例6-6-2】运动声源造成的声波的多普勒效应。这从图形上不大好观察，要转化成声音才能很清楚地感受到。【例7-1-4】进行四连杆机构的运动分析时，由四根连杆的几何位置方程(非线性方程)解出角位置数组，然后经过求导得出角速度和角加速度。【例8-3-1】可控硅电流波形的绘制。特别注意其中出现了不连续点，对不连

22、续点的绘制需要在该点设置左右边界，才能得出正确的图形。【例8-4-2】给出了平面上参数不同的曲线族综合绘制的方法，得出了著名的阻抗圆图。【例9-1-1】典型连续信号波形的绘制。注意要用stairs函数来绘制不连续的阶跃波形。【例9-2-1】典型离散信号波形的绘制。注意要用stem函数来绘制脉冲波形。本节习题本节习题1.(a)在同一张图中画出函数f(x)=sinx+cos(x/2)和g(x)=lnx+x1/3。(b)求函数f(x)的周期。(c)用代数方法证明(b)中的结果。(d)在同一张图中画出h(x)=f(x)-g(x)。2.设x=a cost，y=b sint，0t2。(a)利用参数方程作出

23、它的曲线。(b)消去t，将此方程化为直角坐标中的形式，作出它的图形。3.写出程序，要求能一次画出y=bx(b=0、1、2、3，-x)在x-y坐标中的曲线族，并对所得的图形做出解释。4.对以下题目，用图形来估计极限。(a)(b)(c)2x16xlim42x2231x1)(x3x5xxlim47x9xlim223x5.(a)画出h(x)=x2 cos(1/x)，并估计，根据需要可以把原点附近的图形放大。(b)再画出k(x)=cos(1/x)，比较h和k在原点附近有什么相同点和不同点。6.对下列函数，求和 。(a)(b)(c)f(x)=x2e-xh(x)lim0 xf(x)limxf(x)limxx

24、|x|lnf(x)x1xsinf(x)7.在-4，4区间上作出函数 的图形，研究 和。8.观察函数的图形并求极限。9.通过画图找到，并用代数方法证明之。xx9xxf(x)33f(x)limxf(x)lim0 xx12x1f(x)x-12x1lim1xf(x)limxxxxxf(x)2210.画出下列函数的曲线并判断其极限。(a)(b)(c)(d)11.计算下列函数在给定的各个区间中的平均变化率。(a)f(x)=ex，(i)x-2，0，(ii)x1，3(b)f(x)=lnx，(i)x1，4，(ii)x100，103(c)f(x)=cotx,(i)x/4，3/4,(ii)x/6，/2lnxcotx

25、ln lim0 xx0 x(x)limsinxxxcosx-sinxlim20 xlnxxlim20 x12.参照例5-1-8的几点要求，对以下的曲线的给定点画出割线和切线。(a)f(x)=x+sin2x，x0=/2(b)f(x)=cosx+4sin2x，x0=13.在区间0 x2中画出曲线 ，并在同一幅图上画出曲线 ，先取h=1、0.5、0.1，再取h=-1、-0.5、-0.1，讨论其结果。)2/(1xy hxhxf(x)14.对下列函数:(a)f(x)=x3+x2-x，x0=1(b)f(x)=x1/3+x2/3，x0=1(c)f(x)=sin2x，x0=/2(d)f(x)=x2cosx，x

26、0=/4按要求进行(i)画出y=f(x)以看出函数的全局特点;(ii)在任意点x定义它用任意步长h求得的差商;(iii)令步长h趋于零，判断得到差商的何种公式;(iv)把x取为x0，画出函数曲线及该点的切线;(v)使(iii)中所得公式中的x上下变动，观察它对导数的影响;(vi)画出(iii)中所得公式的曲线，说明它的正、负、零区域意味着什么。5.2 解析几何和多变量分析解析几何和多变量分析5.2.1 平面曲线的绘制平面曲线的绘制【例例5-2-1】绘制极坐标系下的平面曲线=a cos(b+n)并讨论参数a、b、n对曲线的影响。解解:建模绘图的基本方法仍然是先设置自变量数组，按元素群运算的要求列

27、出函数表示式，使得一个表示式能够同时计算出与自变量数目相等的因变量，即可用来绘图。为了便于比较，编一个能同时画出两个图形的程序，采用for循环，读者可从中看到利用循环指数的技巧。由于有两种参数的数据，因变量也有两组数据，因此rho要设成二维数组。MATLAB程序 theta=0:0.1:2*pi;%产生极角向量，这条语句的缺点是3.1到2*pi之间不能被覆盖，改进的方法是用 theta=linspace(0，2*pi，N)，它把0到2*pi之间均分为N份 for i=1:2 a(i)=input(a=);b(i)=input(b=);n(i)=input(n=)rho(i，:)=a(i)*co

28、s(b(i)+n(i)*theta);%极坐标方程 subplot(1，2，i)，polar(theta，rho(i，:);%极坐标系绘图 end 程序运行结果运行并输入不同参数的结果如图5-12所示。a=2;b=pi/4;n=2 (四叶玫瑰线，图5-12(a)a=2;b=0;n=3 (三叶玫瑰线，图5-12(b)图 5-12 例5-2-1中的曲线【例例5-2-2】根据开普勒定律，当忽略其他天体引力时，在一个大天体引力下的小天体应该取椭圆、抛物线或双曲线轨道。它相对于大天体的极坐标位置应满足下列方程:=+cos其中，为偏心率，为常数，为向径，为极角，所取的轨道形状由偏心率决定。试画出天体轨道，

29、并讨论参数、对轨道形状的影响。解解:建模移项，将方程变换为给定和，输入一个数组，就可以得出同样长度的数组，为绘图提供了完整的数据。cos1 MATLAB程序 theta=0:0.1:2*pi;%产生自变量数组 beta=input(beta=);(本书取为1)for i=1:3 epsilon=input(epsilon=);(本书依次取为0.5，1，1.5)rho=beta./(1-epsilon.*cos(theta);%求因变量数组 subplot(1，3，i)，polar(theta，rho);%极坐标系绘图 end 程序运行结果程序运行结果如图5-13所示。其中第三个图比例尺太大，看

30、不出原点附近的情况，因此键入改变图形比例尺的语句:subplot(1，3，3)，axis(-2，2，-2，2)改变将成比例地改变轨道的直径，这是显而易见的。在天体运动中，原点位于椭圆轨道的焦点上，其方程用极坐标来表示比较方便。在平面解析几何中，以原点作为中心的椭圆方程，通常表示为直角坐标的形式:1byax2222图 5-13 不同偏心率时所得的行星轨道(a)椭圆(1);(b)抛物线(=1);(c)双曲线(1)如果在直角坐标中绘图，按照常规，就要把y表示为x的显函数，即这里正负号都要考虑，否则画出的椭圆就少了一半。在给自变量x赋值的时候，x不宜太大，略大于a即可。不要让出现虚数的y的无效数据过多

31、。绘图程序如下：x=linspace(-pi，pi，1001);%自变量数组 y1=2*sqrt(1-x.2/32);%上半部因变量计算 y2=-2*sqrt(1-x.2/32);%下半部因变量计算 plot(x，y1;y2)，grid on%两半图都画 axis equal%使x，y轴同比例22/ax1by最后一条语句是为了保持椭圆形状的纵横轴比例与给定的数据相吻合。得到的曲线如图5-14所示，图中已经对图形线型进行了编辑，把y2代表的负半周椭圆用点画线表示。用参数方程的形式作图也许更为方便，因为它不必考虑开方的正负号和出现的虚数。椭圆的参数方程形式如下:x=a cos，y=b sin因此自

32、变量可以设为，其范围为02。有 theta=linspace(-pi，pi，1001);plot(3*cos(theta)，2*sin(theta)，grid on，axis equal得出的图形是一样的。图 5-14 直角坐标系中的椭圆绘制5.2.2 二次曲面的画法二次曲面的画法【例例5-2-3】二次曲面的方程如下要求讨论参数a、b、c对曲面形状的影响，并画出其图形。解解:建模此题数学模型很清楚，关键在于如何作出三维曲面图形。首先，自变量x和y不再是一维数组，而要设置成二维的矩阵(网格);其次，从上题中知道，在按给定的x、y求z时有开方运算，正负结果都要考虑在内，即按此式计算，一是有正负两个

33、解，这在例5-2-2中已经谈到;二是在x、y取某些值时，z会出现虚数，而在绘制三维曲面的mesh(x，y，z)或surf(x，y，z)函数中，x、y、z只允许为实矩阵，为了使虚数不出现在绘图中，把z的实部z1=real(z)作为三维绘图的因变量。1czbyax2222222222/by/ax1bz MATLAB程序 a=input(a=);b=input(b=);%输入参数，N为网格线的数目 c=input(c=);N=input(N=);xgrid=linspace(-abs(a)，abs(a)，N);%建立x网格坐标 ygrid=linspace(-abs(b)，abs(b)，N);%建立

34、y网格坐标 x，y=meshgrid(xgrid，ygrid);%自变量x，y矩阵 z=c*sqrt(1-y.*y/b/b-x.*x/a/a);%因变量矩阵 z1=real(z);%取z的实部z1(去掉虚数)surf(x，y，z1)，hold on，surf(x，y，-z1);%画正负两半空间曲面标注语句略去。程序运行结果运行这个程序，系统就会提示用户依次输入a、b、c和N，然后给出相应的曲面。将程序稍作修改，使它可以把上述程序循环运行三次，输入的变量第一次为a，b，c，N=5，4，3，20，得到的是椭球;第二次为a，b，c，N=5i，4，3，15，得到的是鞍面;第三次为a，b，c，N=5i，

35、4i，3，10，得到的是双曲面。把画出的图形分别画在三个子图中，其结果如图5-15所示，要特别注意a、b取虚数时对图形的影响，并作讨论。图中的椭球在过x-y平面处有许多“毛刺”，这是由x和y离散取值，不能对准真正的z的过零点造成的。可以用增加网格密度的方法来改进，更好的办法是用参数方程来画椭球曲面。要看出曲面的真实形状，可键入axis equal，使三维尺寸按同样的比例显示。图 5-15 a、b、c、N参数取不同值时所得的不同形状的二次曲面(a)a,b,c,N=5,4,3,20;(b)a,b,c,N=5i,4,3,15;(c)a,b,c,N=5i,4i,3,10【例例5-2-4】与例5-2-3

36、 中的直角坐标方程等价的椭球参数方程为:x=a cosu sinvy=b sinu sinvz=c cosv其中u的取值范围为02，其几何意义是球面上任何点向径的投影在x-y平面上的方位角,即相当于例5-2-2中的。v的取值范围是-+，其几何意义是球面上任何点向径的俯仰角。这个参数方程避开了上下两个半椭球不同方程的衔接问题，使它成为连续平滑过渡的同一个方程，不仅简化了程序，而且画出的图形也漂亮得多。MATLAB程序设a=3，b=2，c=4，则程序为:a=3，b=2，c=4，u=linspace(0，2*pi，20);%将u设为列向量 v=linspace(-pi，pi，20);%将v设为行向量

37、 x=a*cos(u)*sin(v);%x为length(u)length(v)矩阵 y=b*sin(u)*sin(v);%y为与x同阶的矩阵 z=c*ones(size(u)*cos(v);%z也为与x同阶的矩阵 mesh(z，x，y)axis equal 程序运行结果执行这个程序，得到如图5-16所示的曲面，虽然我们对u、v的分割取的是20份，并不是特别细，但是图形还是非常光滑的。在这个程序中有几个问题需要解释清楚:(1)二维参数作为自变量必须构成二维向量，也就是要能构成一个矩阵，而不能是一个向量。因此要由一个自变量的列向量乘以另一个自变量的行向量来构成自变量矩阵。在本程序中，我们将u设为

38、列向量，v设为行向量，x和y的表示式中本来就有u和v的函数的乘积，它们就自然构成了u向量长度和v向量长度相乘阶数的矩阵。图 5-16 用参数方程画出的椭球(2)对于z，它的表达式z=c cosv中只有v，没有u。如果程序中简单地用z=c*cos(v),那得出的z将是与v同长度的行向量,而不是矩阵，与x及y中的元素数目完全不同，无法放在一起画曲面。其实z的表达式中没有u，表示它与u无关，所以用一个阶数与u相同的全么向量ones(size(u)去左乘它，就可以使z的阶数与x及y相同，且不影响z的值。有了这个程序作基础，其他类似的参数曲面就很容易画了。例如方程为/2u2ezusinvyucosvx的

39、曲面就可以用下列程序画出:u=0:0.03:3;v=linspace(0，2*pi);x=u*cos(v);y=u*sin(v);z=exp(-u.*u*ones(size(v)/2);mesh(x，y，z)其所得曲面如图5-17所示。图 5-17 用参数方程画出的帽型曲面5.2.3 空间两曲面的交线空间两曲面的交线【例例 5-2-5】给出空间两曲面方程z1=f1(x，y)及z2=f2(x,y)，列出绘制两曲面及其交线的MATLAB程序。解解:建模两空间曲面方程联立起来，就形成一个空间曲线的方程，这个曲线能同时满足两个曲面的方程，因而也就是这两个空间曲面的交线。显示这两个曲面不难，用两次mes

40、h语句,之间用一条hold on语句即可。要显示其交线，必须先找到各个交点。因为数值计算得到的是离散点，难于找到两曲面上的完全重合的点，本程序采用了设置门限的方法，只要同样网格点处两曲面的z值之差小于设定的门限，就认为它是交点，门限值要试验几次才能定得合适。MATLAB程序%本程序给出两个空间曲面的交线(当然是空间曲线)，给出不同的z1，z2方程%即可画出不同的空间曲面和曲线。x，y=meshgrid(-2:.1:2);%确定计算和绘图的定义域网格 z1=x.*x-2*y.*y;%第一个曲面方程 z2=2*x-3*y;%第二个曲面方程(平面)mesh(x，y，z1);hold on;mesh(

41、x，y，z2);%分别画出两个曲面 r0=abs(z1-z2)=.1;%求两曲面z坐标之差小于0.1的网格矩阵 zz=r0.*z1;yy=r0.*y;xx=r0.*x;%求这些网格上的坐标值，即交线坐标值 plot3(xx(r0=0)，yy(r0=0)，zz(r0=0)，*);%画出这些点不用彩色而用灰度表示 colormap(gray)，hold off 程序运行结果执行此程序得出的曲面如图5-18所示。如果要改变曲面方程，可以在程序中改动第二行和第三行。这还是不好，因为程序还不是通用的。应该使程序运行时向用户提问，允许用户输入任何曲面的方程。此时就要用到字符串功能和eval命令。s1=in

42、put(输入方程语句，s);将原来的z1方程语句改为 z1=eval(s1);图 5-18 曲面与平面的交线这样就可以得出绘制两个任意曲面的交线的程序exn525a。此外，最好让用户能给出定义域和间隔，这比较简单，只要把程序的第一句改为:x，y=meshgrid(xmin:dx:xmax，ymin:dy:ymax);并使xmin、dx、xmax、ymin、dy、ymax都可由人机界面提问并由键盘输入。求交点的门限最好也能人为设定，但这又增加了运行时的麻烦。所以编程时要做出一个折中的选择，即既要有一定的灵活性又要恰到好处。【例例5-2-6】绘制由水平截面与方程z=x2-2y2构成的马鞍面形成的交

43、线，并讨论等高线和方向导数(梯度)的意义。解解:建模对例5-2-5的程序作如下修改:定义域网格改为x，y=meshgrid(-10:.2:10);第一个曲面方程改为z1=(x.2-2*y.2)+eps;第二个曲面(平面)方程改为与z轴正交的水平面z2=a，为了画z2平面，应使z2与x、y有同样的维数，故写成z2=a*ones(size(x)，a可由用户输入。水平平面与曲面的交线还有一个重要的名称，就是等高线。MATLAB三维绘图库中有绘制等高线的命令contour和contour3，前者是把等高线画在x-y平面上，后者则是把等高线画在一定高度的平面上，使之成为立体的，与所在曲面对应。本题中用s

44、ubplot命令把曲面和交线分别画在两幅图上，又把相应的立体化的等高线画在第三幅图上，以便于比较。MATLAB程序 clf，clear x，y=meshgrid(-10:.2:10);%确定计算和绘图的定义域网格 z1=(x.2-2*y.2)+eps;%第一个曲面方程 a=input(a=(-50a50);z2=a*ones(size(x);%第二个曲面方程(平面)subplot(1,3,1),mesh(x,y,z1);hold on;mesh(x，y，z2);%分别画出两个曲面 v=-10 10 -10 10 -100 100;axis(v)，grid%确定第一个分图的坐标系 colorma

45、p(gray)，hold off，%取消彩色改为灰度 r0=abs(z1-z2)0时为极值点;dis0时，若fxxm0，则fm为极小值;若fxxmg(x)的区域，也有f(x)g(x)的区域。若要求两曲线所围的总面积(不管正负)，则可键入 s=trapz(abs(f-g)*dx，在dx=0.001时，得到s=6.47743996919702xcos2)(x2ef(x)图 5-27 被积函数的图形若要求两曲线所围的f(x)g(x)的正面积，则需要一定的技巧，也可有多种方法。(1)方法一 先求出两个交点，根据交点来规定适当的积分上下限。也可只求一个交点，求出左半部分的面积，再乘以2。为了求出交点，可

46、以调用fzero函数，它的第一个变元是函数的表达式，第二个变元是在交点附近的x的初始猜测点。于是键入:x1=fzero(exp(-(x-2).2.*cos(pi*x)-4*cos(x-2)，1)求得 x1=1.06258073077179然后可把积分限按此设定为0 x1，求出半边的积分结果再乘以2，程序如下:x=0 dx x1;f=exp(-(x-2).2.*cos(pi*x);g=4*cos(x-2);s1=2*trapz(abs(f-g)*dx设定dx=0.001时，得到s1=2.30330486000857(2)方法二 设法调用MATLAB中已有的求面积函数quad。这里的关键是要建立一

47、个适当的函数文件，把f(x)-g(x)0的部分取出来。令此函数文件名为exn532f.m，考虑到函数文件必须满足函数群运算的法则，即输入一个x数组，能得出一个同样长度的差值数组，并要使其负差值取零。流程控制语句肯定是不行的，因为它是按单个元素的正负号来决定取舍的，不符合元素群运算规则。此处我们利用了逻辑算式e10，这个函数在e10处取值为1，在e10与e1作元素群乘法，正的e1将全部保留，而负的e1就全部为零了。因此编出的子程序如下:function e=exn532f(x)e1=exp(-(x-2).2.*cos(pi*x)-4*cos(x-2);e=(e1=0).*e1;这个文件应存在MA

48、TLAB的搜索路径下。于是可键入求此积分的主程序语句 s2=quad(exn532f，0，4)得到的结果为 s2=2.30330244952618可以看出，这时求数字积分的关键在于正确地编写待积分函数的子程序。5.3.3 求曲线长度求曲线长度【例例5-3-3】设曲线方程及定义域为 ，-1x1，用计算机做如下工作:(a)按给定区间画出曲线，再按n=2、4、8分割区间并画出割线。(b)求这些线段长度之和，作为弧长的近似值。(c)用积分来估算弧长，并与用割线计算的结果比较，进行讨论。解解:建模先按分区间算割线长度的方法编程，然后令分段数不断增加求得其精密的结果，最后可以与解析结果进行比较。因此编程应

49、该具有普遍性，能由用户设定段数，并在任何段数下算出结果。2x1f(x)MATLAB程序 n=input(分段数目n=)，%输入分段数目 x=linspace(-1，1，n+1);%设定x向量 y=sqrt(1-x.2);%求y向量 plot(x，y)%绘图 hold on%保持 Dx=diff(x);%求各段割线的x方向长度%注意x向量长度为n+1，Dx是相邻两个x元素的差，其元素数为n Dy=diff(y);%Dy是相邻两个y元素的差，其元素数也为n Ln=sqrt(Dx.2+Dy.2);%求各割线长度 L=sum(Ln)%求n段割线的总长度 程序运行结果程序运行后得到图5-28，对于不同的

50、n，其数值结果如下：n=2，L=2.82842712474619 n=4，L=3.03527618041008 n=8，L=3.10449606777484 n=1000，L=3.14156635621648图 5-28 用割线长度逼近弧长我们已经可以大致猜测出它将趋向于，精确的极限值可用符号数学导出。这个程序其实有相当的通用性，不同的被积函数，只要改变其中的一条函数赋值语句，并相应地改变自变量的赋值范围就行了。用符号数学计算此定积分的公式应该是故其相应的程序为 x=syms x,y=sqrt(1-x2)，L=int(sqrt(1+diff(y)2)，-1，1)得出的结果为 L=piba2dx