MATLAB在高等数学中的应用举例（课件PPT） 222页

'第二篇数学篇MATLAB在高等数学中的应用举例 数值解还是解析解？在“语言篇”介绍了MATLAB数学软件基本语法的基础上，本篇讨论如何用它来解决大学工科所学过的数学问题。本篇的第1～4节依次讨论微积分中的极限与导数、解析几何、数列和级数以及数字积分问题；第5节讨论线性代数，第6节讨论概率论与数理统计。主要通过一些例题说明如何灵活使用MATLAB的各种函数来解题。MATLAB是从数值运算发展起来的，以后发展到符号数学（即公式推导）。数学界比较喜欢用Mathematica软件，我们以讨论数值计算方法为主，并考虑到与后续课程的衔接问题，所以选择MATLAB软件。对于工科而言，各门后续课程和未来的工程实践中遇到最大量的将是数值计算问题。计算机首先是计算的工具。计算机的计算过程和方法都是从计算器升级而来的，学生可以理解和接受其每一步，甚至自己都可以编出相应的程序，这是数值计算的一个长处。相比起来，计算机推理过程对于工科学生而言，往往只能知其然而不知其所以然。 数值解还是解析解？(续)其次，用推理方法只能解决很少一部分有解析解的数学命题。比如许多函数是无法求不定积分的，而它们的数值积分却都可以求得。因此优先让学生掌握数值方法等于给他们教会了具有更普遍适用的科学计算方法。对于他们今后的工程生涯将有更广泛的用处。当然，用符号数学求积分，或求微分方程的解析解有时也非常有用，这相当于查积分表或数学手册的电子化和智能化。当一个例题可以同时用数值方法和符号推理方法解决时，我们采取的原则是优先讲数值方法，当同时介绍两种方法时，书中将把符号推理方法的程序用文字框框出来，便于读者跳过不学。因为通常这两种方法在编程上很相似，但又有一些关键性的差别，初学者很容易混淆。对于这类读者，还是先掌握数值方法为好。 5.1函数、极限和导数一．单变量函数值的计算和绘图【例5-1-1】单变量函数的计算和绘图,设要求以0.01s为间隔，求出y的151个点，并求出其导数的值和曲线。解:◆建模:可以采取下列两种方法来做,1。直接用文本文件中编程的方法;2。编成函数文件，由主程序调用的方法;求导数采用diff函数对数组y作运算的方法. 主程序方法主程序exn511a为：t=[0:.01:1.5];%设定自变量数组tw=4*sqrt(3);%固定频率y=w/8*exp(-4*t).*sin(w*t+pi/3);%用数组运算式subplot(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%绘制导数曲线 用函数程序方法此方法,要用两个程序文件主程序exn511b为：dt=0.01;t=[0:dt:1.5];w=4*sqrt(3);y=exn511bf(t,w);Dy=diff(y)/dt;%绘图和加标注的程序略去另要建立一个函数文件exn511bf.m,其内容为：functionxvalues=ex511bf(tvalues,w)%此函数文件,应该能用元素群运算。xvalues=w*exp(-4*tvalues).*sin(w*tvalues+pi/3); 程序运行结果运行这两种程序都得到图5-1-1的曲线。为了节省篇幅，我们没让显示y的数据。以后的各例中还将删略绘图时的标注语句。从本例看，第二种方法似乎更麻烦一些，但它具备模块化的特点。当程序中要反复多次调用此函数，而且输入不同自变量时，利用函数文件可大大简化编程。我们应该掌握这种方法。两次应用diff函数或用diff(y,2)可以求y的二次导数。改变背景色为白色的命令set(gcf,"color","w") 参变方程的计算和绘图【例5-1-2】摆线的绘制当圆轮在平面上滚动时，轮上任一点所画出的轨迹称为摆线。如果这一点不在圆周上而在圆内，则生成内摆线；如果该点在圆外，即离圆心距离大于半径，则生成外摆线。 摆线绘制的程序◆建模：其普遍方程可表为：x=rt-Rsinty=Rcost◆MATLAB程序exn512：t=0:0.1:10;%设定参数数组r=input("r="),R=input("R=")%输入常数x=r*t-R*sin(t);y=r-R*cos(t);%计算x,yplot(x,y),axis("equal")%绘图 摆线绘制程序的结果设r=1,令R=r,R=0.7及R=1.5时得到的摆线、内摆线和外摆线都绘于图5-1-3中。为了显示摆线的正确形状，x,y坐标保持等比例是很重要的，因此程序中要加axis(‘equal’)语句。 三曲线族的绘制【例5-1-3】三次曲线的方程为,试探讨参数a和c对其图形的影响.图5-1-3c和a取不同值时y=ax3+cx的曲线族解:◆方法因为函数比较简单,可以直接写入绘图语句中,用循环语句来改变参数.注意坐标的设定方法,以得到适于观察的图形。给出的程序不是唯一的,例如也可用fplot函数等, 画曲线族的程序x=-2:0.1:2;%给定x数组,subplot(1,2,1)%分两个画面绘图forc=-3:3plot(x,x.^3+c*x),holdon,%a=1,取不同的cend,gridonaxis("equal"),axis([-22-33]),%x,y坐标等比例subplot(1,2,2),fora=-3:3plot(x,a*x.^3+x),%c=1,取不同的aholdon,end,gridonaxis("equal"),axis([-22-33]) 例5-1-3程序运行的结果 四．极限判别【例5-1-4】极限的定义和判别用MATLAB来表达推理过程是比较困难的，因为它必须与实际的数值联系起来。比如无法用无穷小和高阶无穷小的概念，只能用e-10,e-20等数值。不过极限的定义恰恰用了δ和ε这些数的概念，因此不难用程序表达。解:◆建模用函数极限的定义对于函数y=f(x),当任意给定一个正数ε时,有一个对应的正数δ存在,使得当时，其中为绝对值符号，则A就是f(x)在时的极限,如果找不到这样的δ,A就不是它的极限。 判断极限条件的程序%任意给出εwhileflag==1epsilon=input("任给一个小的数ε=")whileabs(A-eval(fxc))>epsilondelta=delta/2;x=xc-delta;%找δ%找不到δ，跳出内循环ifabs(delta)0.0001f=x0.^3+10*x0.^2-2*sin(x0)-50;%求f(x0)g=3*x0^2+20*x0-2*cos(x0);%求导数g(x0)xx=x0-f/g;%切线法求解公式e=abs(xx-x0);x0=xx%把新值赋予x0plot(x0,0,"*"),pause(1)%画出新点的位置end 利用MATLAB库函数的解法运行此程序，设x0分别为-10,-2,2,得出的解分别为-9.4384,-2.5648,2.0707,并可从图上看到解逐渐接近精确解的过程。本程序中用了导数的解析形式，这不是一个值得推荐的程序,好的程序应该用数值方法求导数。读者可自行改进这个程序。实际上可直接用fzero求f(x)的根，方法是把f(x)定义为一个子程序函数ex517f.m，其内容为functiony=ex517f(x)y=x.^3+10*x.^2-2*sin(x)-50;执行exn517b，它先要求用户输入x0的值，然后执行x=fzero(‘ex518f’,x0)，即可得出靠x0较近的数值解x。 割线和切线的绘制【例5-1-8】画出下列曲线在x0处的割线和切线：:设，做以下的工作：（a）在区间内画出；（b）保持不动，求差商：作为h的函数。（c）求h→0时q(h)的极限；（d）对于h=3,2,1，定义并画出此割线。（e）在三个点上分别画出曲线的切线。 割线和切线的绘制(续)解：原理：先在指定的定义域内画出函数的图形，然后在附近计算不同h时的差商及其在h→0时的极限。再按割线公式画出割线图形。◆程序exn518为x=linspace(-1/2,3,100);plot(x,x.^3+2*x),gridon,holdon%画曲线%取不同步长求斜率，画出多根割线。h=[1,0.1,0.001];%取不同步长fork=1:3x1=[1,1+h(k)];f=x1.^3+2*x1;q=diff(f)/diff(x1),%求差商plot(x,f(1)+q*(x-x1(1)),":"),%画割线end 割线和切线的绘制(续)%在多个点上，取h=0.1求斜率，画出割线。fork=0:2x1=[k,k+0.1];f=x1.^3+2*x1;%取不同位置x1q=diff(f)/diff(x1),%求差商subplot(1,2,2),plot(x,x.^3+2*x),gridon,holdonplot(x,f(1)+q*(x-x1(1)),":r"),%画割线end 求极值的方法【例5-1-9】求函数位于区间[0,π]内的近似极值。采用数值方法来解这个问题，画图并求其极值。程序如下。clear,clf,x=linspace(0,pi,10000);y=2*sin(2*x).^2+5/2*x.*cos(x/2).^2;plot(x,y),gridon,holdon程序运行后得到右图的曲线。 求极值的方法1—直接法由y的数据直接求极值(exn519a)：求极值大小：>>ymax=max(y)得为：ymax=3.7323求极值处的x坐标，先用find函数找到变量的下标：%找最大值处y的下标>>nm1=find(y==max(y))%求该下标处的x的值>>xm1=x(nm1)得到：nm1=2752,xm1=0.8643。 直接法求局部极值其他两个极值点是局部极值，只能在划分出的局部区域中寻找。先取两个边界点x01=1和x02=2，求两点处的下标值，并求两点之间的极小值：>>n01=find(abs(x-1)<2e-4);%边界处y的下标>>n02=find(abs(x-2)<2e-4);>>ymin=min(y([n01:n02]))%求x=1~2之间ymin%求y的最小值处的下标及x值>>nm2=find(y([n01:n02])==ymin)+n01-1>>xm2=x(nm2 直接法求局部极值的库函数MATLAB已提供了求极小值的函数fminbnd，可以用一条语句求出极小值（参阅4-4-2节）。求极大值时只需将函数的值反号。因此键入下面三条语句就可代替前面的八条语句。[xm1,ym1]=...fminbnd("-(2*sin(2*x).^2+5/2*x.*cos(x/2).^2)",0,1)[xm2,ym2]=...fminbnd("2*sin(2*x).^2+5/2*x.*cos(x/2).^2",1,2)[xm3,ym3]=...fminbnd("-(2*sin(2*x).^2+5/2*x.*cos(x/2).^2)",2,3) 方法2—用Dy=0特性求极值(2)方法二，本方法则利用导数为零性质，但用数值方法求导，其程序如下(exn519b)：>>h=pi/(10000-1);Dy=diff(y)/h;%先求导数>>plot(x(1:end-1),Dy,"-.r")%绘制导数曲线%绘图时注意导数比原函数的长度小1，检测的门限值需要经过试探。因为x的取值是离散的，不但不可能正好取到Dy=0的点，用户也不知道最靠近零的那个点的Dy值有多大。>>kh=input("过零点检测门限取h的kh倍,kh=")>>nm=find(abs(Dy)>ym=y(nm)%求该点函数值， 用Dy=0特性求极值(续)试探的结果，若kh=1,得出的nm只包含两项，意味着只找到两个过零点，有一个没找到，需加大，取kh=8倍，nm又出现了4项，后两项实际上是一个零点。于是有：>>nm=find(abs(Dy)<8*h),ym=y(nm)，xm=x(nm)n=2751516971457146ym=3.73231.94462.95712.9571xm=0.86401.62372.24462.2449最后一项可以去掉，得到所需的三个极值解。 方法3—用符号数学好处是求过零点不要试凑，用fsolve语句，程序exn519c>>y="2*sin(2*x)^2+5/2*x*cos(x/2)^2";>>Dy=diff(y);>>ezplot(y,[0,pi]),holdon,ezplot(Dy,[0,pi]),gridon根据Dy=0的条件解出相应的xm值。因y在x=1,1.5.2.2附近各有一个极值点,要用三条fsolve语句分别求出:>>x1=fsolve("4*sin(4*x)+5/4*cos(x)+5/4-5/4*x*sin(x)",1)>>x2=fsolve("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]),得:ym=3.73231.94462.9571 应用篇中与本节相关的例题①【例6-2-2】，平面上质点运动轨迹的绘制，即知道参数方程x=x(t)及y=y(t)，画出y=f(x)。②【例6-2-3】平面上质点运动轨迹及李萨如图形的绘制。与上题相仿，其特点在于这两个参数方程都是周期性的函数，因此得出的李萨如图形。③【例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】可控硅电流波形的绘制。特别注意其中出现了不连续点，对不连续点的绘制需要在该点设置左右边界，才能得出正确图形。⑦【例8-4-2】给出了平面上参数不同的曲线族综合绘制的方法，得出了著名的阻抗园图⑧【例9-1-1】典型连续信号波形的绘制。注意要用stairs函数来绘制不连续的阶跃波形。⑨【例9-2-1】典型离散信号波形的绘制。注意要用stem函数来绘制脉冲。 5.2节解析几何和多变量函数 5.2节解析几何和多变量分析本节主要讨论用MATLAB解决空间解析几何的问题。重点是利用MATLAB的三维作图功能，对多变量函数的性质进行更为形象的讨论。读者要从中学会曲面的绘制方法，等高线图和方向导数（梯度）的计算和绘制，进一步把它们的几何意义与解析表示式联系起来。 极坐标系中的绘图【例5-2-1】绘制极坐标系下的平面曲线并讨论参数a、b、n的影响。◆建模绘图的基本方法仍然是先设置自变量数组，按元数群运算的要求列出函数表示式，使得一个表示式能够同时计算出与自变量数目相等的因变量，立即可用来绘图。为了便于比较，编一个能分别画出两个图形的程序，采用for循环，读者可从中看到利用循环指数的技巧。由于有两种参数下的数据，因变量也有两组数据，所以rho要设成二维的数组。 极坐标绘图程序exn521theta=0:0.1:2*pi;%产生自变量(极角)向量%或用theta=linspace(0,2*pi,N);fori=1:2a(i)=input("a=（书上值为a(1)=2;a(2)=2)");b(i)=input("b=（书上值为b(1)=pi/4;b(2)=0)");n(i)=input("n=（书上值为n(1)=2;n(2)=3)")rho(i,:)=a(i)*cos(b(i)+n(i)*theta);%因变量方程subplot(1,2,i),polar(theta,rho(i,:));%绘图end 极坐标绘图程序的输出图形◆运行并输入不同参数的结果如图5-1-2。a=2;b=pi/4;n=2(四叶玫瑰线)a=2;b=0;n=3(三叶玫瑰线) 行星轨道的绘制【例5-2-2】根据开普勒定律，当忽略其他天体引力，在一个大天体引力下的小天体应该取椭圆、抛物线或双曲线轨道。它相对于大天体的极坐标位置应满足下列方程：其中ε为偏心率，β为常数，ρ为向径而θ为极角，所取的轨道形状由偏心率决定。试画出天体轨道，讨论参数β，ε对轨道形状的影响。解：◆建模：移项，将方程变换为其中自变量为theta，因变量为rho。beta和epsilon为参数。 行星轨道绘制程序exn522theta=0:0.1:2*pi;%产生自变量向量beta=input("beta=(取为1)");fori=1:3%设置小于1，等于1，大于1的三种偏心率epsilon=input("epsilon=(依次取为0.5,1,1.5)");rho=beta./(1-epsilon.*cos(theta));%求因变量subplot(1,3,i),polar(theta,rho);%绘图end 行星轨道绘制程序运行的结果其中第三个图比例尺太大，看不出原点附近的情况。因此键入改变图形比例尺的语句：subplot(1,3,3),axis([-2,2,-2,2])改变beta将成比例地改变轨道的直径。 直角坐标系中椭圆的绘制直角坐标中的椭圆方程为：绘图时把y表示为x的显函数这里正负号必须都要考虑，否则画出的椭圆就少了一半。设a=3,b=2,程序可编写如下：x=linspace(-pi,pi,1001);%自变量数组y1=2*sqrt(1-x.^2/3^2);%上半部因变量计算y2=-2*sqrt(1-x.^2/3^2);%下半部因变量计算plot(x,[y1;y2]),gridon%两半图都画axisequal%x,y轴同比例 参数法绘制直角坐标中的椭圆用参数方程的形式作图也许更为方便，因为它不必考虑开方的正负号和出现的虚数。椭圆的参数方程形式如下：因此自变量可以设为theta，其范围为0～2π。theta=linspace(-pi,pi,1001);plot(3*cos(theta),2*sin(theta)),%省略了变量a,bgridon,axisequal得出的图形是一样的。 二．二次曲面的画法【例5-2-3】二次曲面的方程如下要求讨论参数a,b,c对其形状的影响,并画出其图形。解：◆原理数学模型很清楚，关键在于如何作出三维曲面图形。首先是自变量x和y不再是一维数组，而要设置成二维的矩阵(网格)。其次从上题中知道，在按给定x,y求z时有开方运算，有正有负，都要考虑在内：按此式计算，一是有正负两个解，这在上例中已经谈到；二是在x,y取某些值时，z会出现虚数。为了使虚数不出现在surf中，把z的实部z1=real(z)作为三维绘图的因变量。 画圆锥曲面的程序exn523%a,b为实数时得椭球，a为虚数时得鞍面，a,b为虚数时得双曲面，a=input("a=");b=input("b=");%输入参数a,b,c,c=input("c=");N=input("N=");%N为网格线的数目xgrid=linspace(-abs(a),abs(a),N);%建立x网格坐标ygrid=linspace(-abs(b),abs(b),N);%建立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),holdon%画正负两半空间曲面surf(x,y,-z1); exn523参数不同时生成的图形 用参数方程生成的椭球【例5-2-4】与上述等价的椭球参数方程为：其中u的取值范围为0～2π，其几何意义是球面上任何点向径的投影在xy平面上的方位角，即相当于上例中的θ。v的取值范围是-π～+π，其几何意义是球面上任何点向径的俯仰角。这个参数方程避开了上下两半个椭球不同方程的衔接问题，使它成为连续平滑过渡的同一个方程，既简化了程序，画出的图形也漂亮得多。 用参数方程画椭球的程序程序exn524a=3,b=2,c=4,u=linspace(0,2*pi,20)";%u设为列向量v=linspace(-pi,pi,20);%v设为行向量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)axisequal 程序exn524生成的图形执行这个程序所得到的图形如右，虽然我们对u,v的分割取的是20份，并不特别细。但是图形却是非常光滑的。 参数方程绘曲面的要点（1）。二维参数作为自变量必须构成二维向量，也就是要能构成一个矩阵，而不能是一个向量。因此要由一个自变量的列向量乘以另一个自变量的行向量来构成自变量矩阵。在本程序中，我们将u设为列向量，v设为行向量，x和y的表示式中本来就有u和v的函数的乘积，它们就自然构成了u向量长度和v向量长度相乘阶数的矩阵。（2）。对于z，它的表达式中只有v，没有u。如果程序中简单地用z=c*cos(v)，那得出的z将是与v同长度的行向量，不是矩阵，与x及y中的元素数目完全不同，无法放在一起画曲面。其实z的表达式中没有u,表示它与u无关，所以用一个阶数与u相同的全么向量ones(size(u))去左乘它，就可以使z的阶数与x及y相同，且不影响z的值。 参数方程绘曲面另一程序%程序exn524a：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)axis([-3,3,-3,3,-1,1]) 三．空间两曲面的交线【例5-2-5】给出空间两曲面方程z1=f1(x,y)及z2=f2(x,y)，列出绘制两曲面及其交线的MATLAB程序。解：◆方法两空间曲面方程联立起来，就形成一个空间曲线的方程，这个曲线能满足两个曲面的方程，因而也就是这两个空间曲面的交线。显示这两个曲面，可用两次mesh语句，之间有一条holdon语句即可。要显示其交线，必须先找到各个交点。因为数值计算得到的是离散点，难于找到两曲面上的完全重合的点，本程序采用了设置门限的方法，只要同样网格点处两曲面的z值之差小于设定门限，就认为它是交点，门限值要试验几次才能定得好 绘制两曲面交线的程序exn525%本程序给出两个空间曲面的交线，方程%给出不同的z1,z2方程，画出不同的空间曲面和曲线。[x,y]=meshgrid(-2:.1:2);%确定定义域网格z1=x.*x-2*y.*y;%第一个曲面方程z2=2*x-3*y;%第二个曲面方程%分别画出两个曲面（保持坐标比例不变）mesh(x,y,z1);hold;mesh(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),’*’);%画出这些点 程序exn525生成的图形 通用的画曲面交线程序exn525a如果改变曲面方程，可以在程序中改动第二和第三行。这究竟不好，程序还不能说是通用的。应该使程序运行时向用户提问，允许用户输入曲面的方程。此时就要用到字符串功能和eval命令。s1=input(‘输入方程语句’,’s’);在原来的z1方程语句处改为z1=eval(s1);这样就可以得出绘制两个任意曲面的交线的程序exn525a。此外，最好让用户能给出定义域和间隔。这比较简单，只要把第一句改为：[x,y]=meshgrid(xmin:dx:xmax,ymin:dy:ymax); 等高线和方向导数（梯度）的绘制【例5-2-6】画出多个水平截面与马鞍面的交线，讨论它们的意义。解:◆方法对上例的程序作如下修改。定义域网格改为[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，前者是把等高线画在xy平面上，后者则是把等高线画在该高度的平面上，使之具有立体的、与所在曲面对应的部位。本题中用subplot命令把曲面和交线分别画在两张图上，又把相应的立体化的等高线画在第三张子图上，以便于比较。， 绘制等高线的程序exn526[x,y]=meshgrid(-10:.2:10);%确定定义域网格z1=(x.^2-2*y.^2)+eps;%第一个曲面方程a=input("a=(-50>x0=roots([3,6,-9])x0=-3.0000,1.0000>>y0=roots([-3,6,0])y0=02知极值点应在A(-3,0),B(-3,2),C(1,0),D(1,2)四处， 多变量函数的极值程序exn526%画出[-4,4][-4,4]的区域内曲面f(x,y)的形状。[x,y]=meshgrid(-4:.1:4);f=x.^3-y.^3+3*x.^2+3*y.^2-9*x,mesh(x,y,f)%计算A,B,C,D四点的函数fm及其各阶偏导数：xm=[-3,-3,1,1];ym=[0,2,0,2];fm=xm.^3-ym.^3+3*xm.^2+3*ym.^2-9*xmfxm=3*xm.^2+6*xm-9,fym=-3*ym.^2+6*xm,fxxm=6*xm+6,fyym=-6*ym+6,fxym=0,dis=fxx.*fyy-(fxy.^2) 例【5.2.7】的空间曲面 多变量函数的极值(续3)最后一项dis是极点性质的判别式，dis>0时为极值点。Dis<0时为鞍点，dis=0时不定。在dis>0时，若fxxm>0，fm为极小值，若fxxm<0，fm为极大值。运行此程序exn527得到fm=[27,31,-5,-1]，fxxm=[-12,-12,12,12]，fyym=[6,-6,6,-6]，fxym=0。判别式dis==[-72,72,72,-72]，说明只有B,C点是极值点。B点两个方向的二阶导数均为负，是极大值31；C点两个方向的二阶导数均为正，是极小值-1。A,D两点则是鞍点。 方法二：等高线法：程序exn527a[x,y]=meshgrid([-5:0.5:3],[-3:0.5:5]);z=x.^3-y.^3+3*x.^2+3*y.^2-9*x;%画出等高线图和彩色标尺图contour(x,y,z,20),colorbar,holdon,gridonplot([-3,1,1,-3],[0,0,2,2],"o")%以0.2步长求z1的梯度的x,y分量[px,py]=gradient(z,.2,.2);%画出梯度向量图quiver(x,y,px,py) 例5.2.7的等高线和梯度图 程序exn527a执行的结果讨论◆程序运行后得到图5-2-7b。在彩色屏幕上，用颜色表示等高线值的大小，旁边的彩色标尺标示了不同的颜色所对应的z值。在黑白印刷的书上，只好藉助于梯度向量图。在极值点，梯度为零，在极值点的附近，等高线应围绕它，呈闭合形式而梯度箭头应都指向该点（极大值），或者背离该点（极小值）。如果某点梯度为零，但等高线不是环绕而是通过它，而既有部分梯度向量指向它，又有部分梯度向量背离它，那它就不是极值点而是鞍点。由这些准则，很容易看出B处是极大值位置而C是极小值位置，而A,D两点则是鞍点。 五．柱面的绘制通常要画的柱面是平行于某个坐标轴的，这时方程中将不出现这个变量，例如在三维空间中方程表示的是平面曲线沿z向无限延伸的柱面。在方程中没有在z，此时三维曲面的绘制需要一点技巧，可从下面的例中看到。【例5-2-8】画出下列方程所表示的三维空间中的柱面：(a)，(b)， 小题(a)的自变量矩阵设置(a)X是一个自变量矩阵，而Z则是因变量，所以另一个自变量矩阵必定是Y。自变量平面是xy平面,其中x是真正的自变量而y轴是柱面方向。设定列向量u为x轴方向，行向量v为y轴方向。X矩阵由列向量u与一个长度同v的全么行向量相乘构成，Y矩阵由长度与u向量相同的全么列向量与v行向量相乘构成，生成的X,Y取值用图5-2-9a给出，要注意的是此处的横坐标恰好对应于X,Y中的列向量，图形和矩阵排列的位置恰好转了90度。而Z则由x与y之间的函数关系确定。 解(a)的程序语句如下:u=linspace(-5,5,10)";v=linspace(-5,10,10);X=u*ones(size(v));Y=ones(size(u))*v;Z=2-X.^2;subplot(1,2,1),mesh(X,Y,Z)此程序运行的结果得到图5-2-8b中的左图。它是一个沿y轴的柱形曲面。 小题(b)的自变量矩阵设置X是一个自变量矩阵，而Y则是因变量，所以另一个自变量矩阵必定是Z。自变量平面是xz平面,而z轴是柱面方向。X矩阵由列向量u与一个长度同v的全么行向量相乘构成，Z矩阵由长度与u向量相同的全么列向量与v行向量相乘构成，而Y则由x与y之间的函数关系确定。因为函数关系的开方包含正负两个分量，要分别用Y1,Y2表示，并分别画图。 解(b)的程序语句如下:X1=u*ones(size(v));Y1=2*sqrt(1+X1.^2);Y2=-2*sqrt(1+X1.^2);Z1=ones(size(u))*v;subplot(1,2,2),mesh(X1,Y1,Z1),holdon,mesh(X1,Y2,Z1)它的运行的结果得到图5-2-8b中的右图。它是一个沿z轴的柱面。 【应用篇中与本节相关的例题】①【例6-4-3】电磁场的分布是一个空间问题，用手工进行计算和绘制比较困难。本题说明了如何计算和绘制空间电位分布的立体图，说明等电位线的计算及绘制，电场向量是等电位线的梯度，所以也可得出电场的分布图。②【例6-5-1】，本例讨论由电流产生的空间磁场计算问题，除了计算公式不同外，磁场与电场问题相仿，只不过其重点放在磁场的空间分布方面。③【例6-5-2】本例继续讨论由双线圈通电流产生的空间磁场计算问题，研究等磁场强度的区域。读者可以看到深入探讨磁场分布及均匀性分析的一些工程表述方法。④【例6-7-1】，【例6-7-2】光波的波形叠加造成的干涉现象。光的干涉也是一个空间问题。为了简化，这些例题还都限于两维空间。而把第三维——光的强度用灰度来表示，从而鲜明地看出干涉的效果。 【应用篇中与本节相关的例题】⑤【例7-1-1】和【例6-4-3】，【例6-5-1】都涉及到空间向量的分析计算，从这些例题可以知道，如何用MATLAB来进行空间向量分析。⑥【例7-1-4】，给出了机械四连杆平面运动中的动画演示。四连杆平面运动的求解是一个非线性问题，求一个点的解就得要靠多次试凑，非常麻烦。要求出许多点就更繁，MATLAB在此处已显示出极大的优越性。现在还要画出它的运动并连贯地以动画显示，又用到了MATLAB强大的绘图功能，所以这一例题确实是手算无法比拟的。⑦【例8-3-2】三相电机旋转磁场的形成过程涉及空间和时间两方面的问题，三个磁极在空间上相差120度，而馈入三个磁极的交流电流在时间上又相差120度，从而形成了旋转磁场。这是很难凭空想象的，本例给出了空间和时间差的综合动画演示，将概念表达得极其清楚，可以说是非常的具有创意的动画演示。 5.3节数列和级数 一．数列的表示方法数列就是自变量为整数时的函数。MATLAB中的元素群运算特别适合于简明地表达数列，可省去其他语言中的循环语句。下面就是一些例子：n=1:6;1./n=1.00000.50000.33330.25000.20000.1667(-1).^n./n=-1.00000.5000-0.33330.2500-0.20000.16671./n./(n+1)=0.50000.16670.08330.05000.03330.0238左端的算式表示这个数列产生方法的“通项”，它必须符合元素群运算的规则，所以要充分注意用点乘、点除和点幂。例如(-1).^n就是产生交项数列符号位的算式，它在n取偶数时为正，而它在n取奇数时为负。在某些情况下，当产生数列的运算中包含数组运算时，就不可避免地要用for循环。 数列用for循环的表示方法比如计算n!（n的阶乘），它应该写成prod(1:n),其中的n就不能是数组，因为prod(1:n)中已用了数组[1:n]。这时必须用：fork=1:6x(k)=1/prod(1:k);end,x得x=1.00000.50000.16670.04170.00830.0014在MATLAB中数列随n增加而变化的趋向很容易由计算其数值并作图的方法来解决。但要求数列在n趋向∞时的极限时往往要藉助于符号数学，可以从下面的实例看出。 【例5-3-1】对下列各题的序列，问：(i)。计算并画出其前25项，判断它是否收敛。若收敛，极限L是多少？(ii)。如果序列收敛，找到数N，使得n>N后的an都有。如果要离极限L小于0.0001，序列该取多长？(1)，(2)，(3)，(4)， 解例5.3.1的程序解：◆只要会写通项的表达式，程序是很简单的。用数值计算方法时，四个题可编在一起如下：程序exn531n=1:25;a1=n.^(1./n);a2=(1+0.5./n).^n;a3=sin(n);a4=n.*sin(1./n);plot(n,a1,n,a2,n,a3,n,a4) 程序exn531的运行结果得到的数列图如右。在计算机屏幕上，四根曲线将用不同的颜色区分和标注，在黑白印刷的书上只好另加字母。可以初步判断，除了a3以外，其他三组数列在n趋向于∞时都趋向于某极限L1,L2,L4。 用符号数学求数列的极限求极限最好用符号数学来解，主要的不同是自变量n应设为符号变量，所有的函数也要重写一次，使它们也成为符号因变量，最好是在程序开始处用clear命令清除掉前面程序在工作空间中生成的同名数值变量。语句如下：clear,symsnL1=limit(n^(1/n),inf)%为了缩短语句，也可写成两句：a1=n^(1/n),L1=limit(n^(1/n),inf)L2=limit((1+0.5./n).^n,inf)L4=limit(n.*sin(1./n),inf)程序运行后，得到L1=1,L2=exp(1/2),L4=1 二．常数项级数无穷数列的累加称为级数，当取其前面若干有限项时，得到的是部分和。将数列a累加形成的新序列可用s=cumsum(a)实现，如果a的长度是n，则s的长度也是n。即每一个s(k)是数组a中前k项的和。注意cumsum与sum命令的区别，若ss=sum(a)，得到的是一个数，是序列s=cumsum(a)中最后一项ss=s(n)。因为它是把a中所有元素加在一起得到的最后结果。MATLAB中同样有符号数学的累加命令，要注意它与数值计算的差别，主要是符号数学没有数组累加成数组的命令，只有求一个求总和的累加命令symcum。 【例5-3-2】设级数（a）,（b），试观察它们的部分和序列变化的趋势，如果是收敛的，计算出它们在n趋向于无穷大时的极限值。解：(1)。用数值方法计算的程序exn532如下：clear,n=input("n=");k=1:n;a1=1./k.^2;s1=cumsum(a1);a2=1./k;s2=cumsum(a2);plot(k,s1,k,s2),gridons1(end),s2(end) 程序exn532的运行结果键入n=20时，得到图形如图，数值结果为：s1(end)=1.59616324391302s2(end)=3.59773965714368我们只能从图形上猜测s1会趋向于一个极限，而s2就难说了。 用符号数学求部分和的程序clear,symsk,ss1_20=symsum(1/k^2,1,20)ss1=symsum(1/k^2,1,inf)ss2=symsum(1/k,1,inf)结果为：ss1_20=1.59616324391302ss1=1/6*pi^2=1.64493406684823ss2=inf 三.函数项级数【例5-3-3】利用幂级数计算指数函数解：◆原理：指数可以展开为幂级数其通项为x^n/prod(1:n)，因此用下列循环相加程序就可计算出这个级数;%输入原始数据，初始化yx=input("x=");n=input("n=");y=1;%将通项循环相加n次，得yfori=1:ny=y+x^i/prod(1:i);end,y分别代入x=1,2,4,-4四个数，取n=10，结果如下表 用级数计算指数的结果exp(x)nx=124-412.000000003.000000005.00000000-3.0000000022.500000005.0000000013.000000005.0000000032.666666676.3333333323.66666667-5.6666666742.708333337.0000000034.333333335.0000000052.716666677.2666666742.86666667-3.5333333362.718055567.3555555648.555555562.1555555672.718253977.3809523851.80634921-1.0952381082.718278777.3873015953.431746030.5301587392.718281537.3887125254.15414462-0.19223986102.718281807.3889947154.443104060.09671958有效数位7420 此计算程序的缺点与问题可以看出这个简单程序虽然原理上正确，但不好用，精度差别很大。问题是：（1）只能用于单个标量x的计算，不能用于x的数组运算；（2）当x为负数时，它成为交项级数，它收敛很慢；（3）此程序要做n2/2次乘法，n很大时，乘法次数太多，计算速度低。（4）若x较大，n就要较大才能达到精度要求。因此n由用户来输入不科学，应该由软件按精度要求来选； 程序改进的方法（一）（1）。考虑到数组和输出显示的程序如下：x=input("x=");n=input("n=");%输入x,n,y=ones(size(x));%初始化yfori=1:ny=y+x.^i/prod(1:i);%循环相加end执行此程序并输入x=[1,2,4,-4]及n=10，可一次得出上表的计算结果。（2）。此时可以利用exp(-x)=1/exp(x)来避免交项级数的计算； 程序改进的方法（二）（3）。设一个中间变量z,它的初始值为z=ones(size(x)),把循环体中的计算语句改成y=y+z;z=x.*z/i这样求得的z就是z=x.^i/i!，于是每个循环只需作一次乘法，计算整个级数只需n次乘法。按这种算法，y的初始值应改为y=zeros(size(x))。（4）。为了按精度选择循环次数，就不该用for循环，而该用while语句，它可以设置循环继续的条件语句。通常可取y+z-y>tol,tol是规定的允许误差。只要相邻两次的y值之差大于tol，循环继续进行，直至小于tol为止。 X太大时exp(x)程序的改进为了使x不致太大，还可以利用关系式exp(x)=（exp(x/k))k，令x1=x/k,k通常取大于而最靠近x的2的幂。（即k=nextpow2(x)）例如x=100,就取k=128,这样保证x1的绝对值小于1，级数收敛得很快。取十项保证有7位有效数。而exp(x1)128可化成x=(...((exp(x1))2)2...)2,即x1的七次自乘。用七次乘法就可完成。这既保证了精度，又提高了速度。 不同阶数泰勒级数的近似程度【例5-3-4】把一个多项式用泰勒级数表示,分析阶数对逼近程度的影响.解:◆原理一个多项式函数可以精确地用泰勒公式展开,但必须取足够高的阶数(等于多项式的次数),否则就会产生误差.在MATLAB中,多项式可以用其系数向量来表示,求值和求导用polyval和polyder命令可参阅本书4.3节。设多项式则它在附近展开的n阶泰勒公式为这个公式是精确的，没有误差项： 泰勒级数展开程序exn534此程序最高只到三阶,如果输入多项式系数向量a的长度不大于4,则其高阶泰勒展开式完全精确,如length(a)大于4,就会有误差,读者可自行试算.并考虑如何编写更完美的程序.a=input("输入多项式系数向量a=[]=");x0=input("展开点的坐标值x0=");[xm]=input("展开坐标区间[xmin,xman]=[]=");x=linspace(xm(1),xm(2));%设定自变量数组y=polyval(a,x);ya=polyval(a,x0);%求y在x0点的值y(x0) 多项式泰勒展开程序(续)Da=polyder(a),Dya=polyval(Da,x0);%求x0点的一阶导数D2a=polyder(Da),D2ya=polyval(D2a,x0);%求x0点的二阶导数D3a=polyder(D2a),D3ya=polyval(D3a,x0);%求x0点的三阶导数yt(1,:)=ya+Dya*(x-x0);%一阶泰勒展开%二阶泰勒展开yt(2,:)=yt(1,:)+D2ya*(x-x0).^2/prod(1:2);%三阶泰勒展开yt(3,:)=yt(2,:)+D3ya*(x-x0).^3/prod(1:3);plot(x,y,’.’,x,yt(1:3,:)),grid%绘图,准确值用点线表示 程序xn534运行结果输入a=[2,-3,4,5];x0=1;xmin=0;xmax=2;得出图5-3-1所示的三根曲线,本来应有四根,但有两根曲线是重合的,因为我们输入的多项式系数向量长度为4.读者可试验输入更长的a来比较其结果. 例5-3-5任意函数的泰勒级数编写演示任意函数展开为各阶泰勒级数的程序,并显示其误差曲线.解:◆原理任意函数的泰勒展开式如下：其中为余项，也就是泰勒级数展开的误差。 泰勒展开程序exn535(1)fxs=input("输入y=f(x)的表达式;%fxs是字符串K=input("输入泰勒级数的展开阶数K=(书上为5)");a=input("展开的位置x0=(书上为0.5)");b=input("展开的区间半宽度b=(书上为2)");x=linspace(a-b,a+b);%构成自变量数组lx=length(x);dx=2*b/(lx-1);%数组x的长度和间距y=eval(fxs);%求出y的准确值%y的准确曲线用点线绘出subplot(1,2,1),plot(x,y,"."),holdon 泰勒展开程序exn535(2)%求出y在a点一阶导数,注意数组求导后长度减一Dy=diff(y)/dx;Dya(1)=Dy(round((lx-1)/2));%一阶泰勒展开,绘图yt(1,:)=y(round(lx/2))+Dya(1)*(x-a);plot(x,yt(1,:))%求出a点2~k阶导数和y的k阶展开fork=2:KDy=diff(y,k)/(dx^k);Dya(k)=Dy(round((lx-k)/2));yt(k,:)=yt(k-1,:)+Dya(k)/prod(1:k)*(x-a).^k;plot(x,yt(k,:)),e(k,:)=y-yt(k,:);%绘图，求出yt的误差end,grid,holdoff,subplot(1,2,2)fork=1:Kplot(x,e(k,:)),holdon,end%绘制误差曲线 程序exn535运行结果输入fxs=cos(x),K=5,a=0.5,b=2所得曲线如下图.读者可改变其坐标系范围以仔细观察最关心的部分。可以看出，泰勒展开的阶数愈高，它与原来函数在愈宽的自变量区间内误差就愈小，即愈逼近于原来的函数。 程序exn535生成的图形 符号数学泰勒展开命令taylor对于本题，其调用格式如下：symsx,f1=taylor(cos(x),8)得到：f1=1-1/2*x^2+1/24*x^4-1/720*x^6第二个变元是展开式的阶数。这样得出的是cos(x)的8阶麦克劳林级数的解析式。如果第二个输入变元a不是整数，它就代表在x=a附近展开的泰勒级数。此时阶数就规定为5，不能由用户指定了。如键入symsx,f1=taylor(cos(x),8.1) 【应用篇中与本节相关的例题】①【例9-1-5】演示了如何用傅立叶级数合成一个周期性的方波，说明傅立叶级数的阶数愈高，其合成波形愈接近于原函数波形。例题并解释了在原函数的间断点附近出现的吉布斯效应。 5.4节数值积分和微分方程数值解 一．数值定积分求面积【例5-4-1】用数值积分法求由，y=0,x=0与x=10围成的图形面积，并讨论步长和积分方法对精度的影响。解:◆原理用矩形法和梯形法分别求数值积分并作比较，步长的变化用循环语句实现。MATLAB中的定积分有专门的函数QUAD,QUADL等实现。为了弄清原理，我们先用直接编程的方法来计算，然后再介绍定积分函数及其调用方法。设x向量的长度取为n，即将积分区间分为n-1段，各段长度为。算出各点的，则矩形法数值积分公式为： 矩形和梯形定积分公式梯形法的公式为：比较两个公式，它们之间的差别只是。在MATLAB中，把向量中各元素叠加的命令是sum。把向量中各元素按梯形法叠加的命令是trapz。梯形法的几何意义是把被积分的函数的各计算点以直线相联，形成许多窄长梯形条，然后叠加，我们把两种算法都编入同一个程序进行比较。 求面积的数值积分程序exn541fordx=[2,1,0.5,0.1]%设不同步长x=0:.1:10;y=-x.*x+115;%取较密的函数样本plot(x,y),holdon%画出被积曲线并保持x1=0:dx:10;y1=-x1.*x1+115;%求取样点上的y1%用矩形（欧拉）法求积分，注意末尾去掉一个点n=length(x1);s=sum(y1(1:n-1))*dx;q=trapz(y1)*dx;%用梯形法求积分stairs(x1,y1),%画出欧拉法的积分区域plot(x1,y1)%画出梯形法的积分区域[dx,s,q],pause(1),holdoff，end 程序exn541运行结果程序运行的结果如下：步长dx矩形法解s梯形法解q29108101865815.5841.25816.25.1821.65816.65用解析法求出的精确解为2450/3=816.6666...。dx=2时矩形法和梯形法的积分面积见图5-4-1.。在曲线的切线斜率为负的情况下，矩形法的积分结果一定偏大，梯形法是由各采样点联线包围的面积，在曲线曲率为负（上凸）时，其积分结果一定偏小，因此精确解在这两者之间。由这结果也能看出，在步长相同时，梯形法的精度比矩形法高。 矩形法数字积分的演示程序rsumsMATLAB中有一个矩形法数字积分的演示程序rsums，可以作一个对比。键入rsums("115-x.^2",0,10)就得到右图。图中表示了被积函数的曲线和被步长分割的小区间，并按各区间中点的函数值构成了各个窄矩形面积。用鼠标拖动图下方的滑尺可以改变步长的值，图的上方显示的是这些矩形面积叠加的结果。 MATLAB内的数值定积分函数在实际工作中，用MATLAB中的定积分求面积的函数quad和quadl可以得到比自编程序更高的精度，因为quad函数用的是辛普生法，即把被积函数用二次曲线逼近的算法，而quadl函数采用了更高阶的逼近方法。它们的调用格式如下：Q=QUADL(FUN,A,B,TOL)其中，FUN是表示被积函数的字符串，A是积分下限，B是积分上限。TOL是规定计算的容差，其默认值为1e-6例如，键入S=quad("-x.*x+115",0,10)得到S=8.166666666666666e+002 二．求两条曲线所围图形的面积【例5-4-2】。设计算区间[0,4]上两曲线所围面积。解：◆原理：先画出图形，>>dx=input("dx=");x=0:dx:4;>>f=exp(-(x-2).^2.*cos(pi*x));>>g=4*cos(x-2);>>plot(x,f,x,g,":r")得到右图。从图上看到，其中既有f(x)>g(x)的区域，也有f(x)>g(x)的区域， 求两条曲线所围图形的面积(1)若要求两曲线所围总面积（不管正负），则可加一条语句>>s=trapz(abs(f-g))*dx,在dx=0.001时，得到s=6.47743996919702若要求两曲线所围的f(x)>g(x)的正面积，则需要一定的技巧.◆方法一。先求出交点x1，再规定积分上下限。>>x1=fzero("exp(-(x-2).^2.*cos(pi*x))-4*cos(x-2)",1)%把积分限设定为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。这里的关键是建立一个函数文件，把e1=f(x)-g(x)>0的部分取出来。利用逻辑算式(e1>0)，它在e1>0处取值为1，在e1<0处则为零。让逻辑函数(e1>0)与e1作元素群乘法，正的e1将全部保留，而负的e1就全部为零。因此编出子程序exn542f.m如下：functione=exn542f(x)e1=exp(-(x-2).^2.*cos(pi*x))-4*cos(x-2);e=(e1>=0).*e1;将它存入工作目录下。于是求此积分的主程序语句为：>>s2=quad("exn542f",0,4)得到的结果为：s2=2.30330244952618 三．求曲线长度【例5-4-3】设曲线方程及定义域为：用计算机做如下工作：(a)按给定区间画出曲线，再按n=2,4,8份分割并画出割线。(b)求这些线段长度之和，作为弧长的近似值。(c)用积分来估算弧长，并与用割线计算的结果比较。解：◆原理：先按分区间算割线长度的方法编程，然后令分段数不断增加求得其精密的结果，最后可以与解析结果进行比较。因此编程应该具有普遍性，能由用户设定段数，并在任何分段数下算出结果。 求曲线长度的程序exn543n=input("分段数目n="),%输入分段数目x=linspace(-1,1,n+1);%设定x向量y=sqrt(1-x.^2);%求y向量plot(x,y),holdon%绘图并保持Dx=diff(x);%求各段割线的x方向长度%x向量长度为n+1，Dx是相邻x元素的差，其元素数为nDy=diff(y);%Dy是相邻两个y元素的差Ln=sqrt(Dx.^2+Dy.^2);%求各割线长度L=sum(Ln)%求n段割线的总长度 程序exn543的运行结果程序运行后得到图5-32，在不同的n下，其数值结果为：n=2,L=2.82842712474619n=4,L=3.03527618041008n=8,L=3.10449606777484n=1000L=3.14156635621648我们已经可以大致猜测出它将趋向于π，精确的极限值可用下列符号数学语句导出。>>x=symsx,y=sqrt(1-x^2),L=int(sqrt(1+diff(y)^2),-1,1)这个程序其实有相当的通用性，不同的被积函数，只要改变其中的一条函数赋值语句，并相应地改变自变量的赋值范围就行了。 四．求旋转体体积【例5-4-4】求曲线与x轴所围成的图形分别绕x轴和y轴旋转所成的旋转体的体积。解：原理：由于旋转对称性，在圆周方向的计算只要乘以圆周长度，不需要积分运算。因此旋转体的体积计算实际上就退化单变量求积分。程序如下：%先画出平面图形>>dx=input("dx=");x=0:dx:pi;>>g=x.*sin(x).^2;plot(x,g) 求筒形旋转体体积(a)。绕y轴体积用薄圆柱筒形体作为微分体积单元，其半径为x，厚度为dx，高度为g(x)，其立体图见图5-34左，此筒形单元的截面积为g*dx，薄环的微体积为：>>dv=2*pi*x*dx*g，旋转体的体积为微分体积单元沿x方向的和，键入：>>v=trapz(2*pi*x.*g*dx)得：v=27.53489480036561 求盘形旋转体体积(b)。绕x轴体积它绕x轴旋转形成一个薄圆盘，其厚度为dx，而半径为g(x)。所以此薄盘体的微体积为：>>dv1=pi*g.^2.*dx，旋转体的体积为微分体积单元沿x方向的和：>>v1=trapz(pi*g.^2*dx)得：v1=9.86294784774497 用符号数学工具箱的程序精确的理论结果可用符号数学工具箱函数求得如下：>>symsx,g=x*sin(x)^2;>>v1t=int(pi*g^2,0,pi),double(v1t)v1t=1/8*pi^4-15/64*pi^2ans=9.86294784774499大多数的定积分并不会有理论的解析结果，所以这样的验证一般是不必要的。 五．多重积分【例5-4-5】计算二重积分积分区域Ω为由x=1,y=x及y=0所围成的闭合区域.解:◆原理先画出积分区域，在任意x处取出沿y向的一个单元条，其宽度为dx，而高度为y=x，所以y是一个数组。其上的被积函数f也是一个数组，沿y向的积分可用trapz(f)完成，得到s1(k)，它是随x而变的。用for循环求出所有的s1(k)。再沿x方向用trapz函数积分。MATLAB的数组运算可以代替一个for循环，所以二重积分只用了一组for语句。 二重积分的MATLAB程序exn545clear,formatcompactfill([0,1,1,0],[0,0,1,0],"y"),hold%画出积分区域fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],"r")%画出单元条dx=input(‘步长dx=‘);dy=dx;x=0:dx:1;lx=length(x);fork=1:lxx1=(k-1)*dx;y1=0:dy:x1;f=x1.^2+y1.^2;s1(k)=trapz(f)*dy;ends=trapz(s1)*dx 用MATLAB函数求二重积分(1)运行的数值结果在步长dx=0.01时为:s=0.3334另一种方法是利用MATLAB中现成的二重积分函数dblquad，其调用格式为：Q=DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX,TOL)其中FUN是x,y的函数，接下来的四个变元是四个积分限，其中前两个对应于x，后两个对应于y，TOL为允许误差（默认值为1.e-16）。这四个积分限只许用常数代入，可见dblquad函数只能用于积分区域为矩形的情况。解决的方法之一是仍用矩形区积分，但把不属于积分区域内的函数置成零，其方法与上题有些类似。 用MATLAB函数求二重积分(2)在图示的积分区域中，对角线左上方的白色区域满足y-x>0，逻辑式(y-x<0)在此区域均等于零,而在灰色区域内为一。将它与被积函数作元素群相乘，就构成了一个新的被积函数，它与原来函数的差别是把灰色积分区外的函数置为零。这样就可以按矩形区域调用dblquad函数了。键入：>>Q=dblquad("(x.^2+y.^2).*(y-x<0)",0,1,0,1)得到Q=0.33332245532028 三重积分的计算【例5-4-6】计算三重积分积分区域Ω为由x=1,y=x,z=xy及三个坐标面所围成的区域.解:◆方法先画出积分区域图5-36,这个区域在xy平面上的投影与图5-35相仿，只是增加了z方向的高度，从而构成了一个三维的实体。先画出积分区域。这个区域在xy平面上的投影上例相仿，只是增加了z方向的高度，从而构成了一个三维的实体。程序exn546a用来画这个立体空间。 x=1,y=x都是沿z向的柱面，本题还用了plot3命令以画出与z轴平行的辅助平面。顶面则是由z=xy构成的二次曲面，用mesh函数容易画出。难点在于此区域有效的xy底面是一个三角形,不易用自变量网格表示。为了解决这个问题，采取了上题中对因变量乘以逻辑式的处理方法。将z乘以(y-x<0)就可使y=x直线左上方所有z值都变成0。画出三维图后可以靠鼠标来拖动三维图形旋转,以得到一个最好的视觉效果。 绘制积分区域的程序exn546a%本程序给出由x=1,y=x,z=xy三个曲面围成的积分区域.[x,y]=meshgrid(0:.05:1);%确定矩形定义域网格z1=x.*y.*(y-x<0);%求z1=xy并构成三角形定义域mesh(x,y,z1);holdon;%画出积分区顶部%以下画出积分区域的几个侧柱面x1=[0:0.02:1];y1=x1;sx1=length(x1);zd=[zeros(1,sx1);x1.*y1];plot3([x1;x1],[y1;y1],zd,"*")line(ones(2,sx1),[y1;y1],[zeros(1,sx1);y1])plot3(ones(2,sx1),[y1;y1],[zeros(1,sx1);y1],"o") 编写三重积分程序的思路(1)在任意点（x,y）处取出沿z向的一个单元条,其底面积为dx*dy,而高度为z=x*y,这一个细柱体上从z=0到z=x*y间的所有各点度属于积分的区域，把它表为z向的一个数组。因为此处(x,y)固定，其上的被积函数f=f(z)是随z而变的一个数组，沿z向的积分可用trapz(f(z))完成，得到s1。(2)s1是随x,y而变的,先固定x,用for循环求出沿y向所有的s1(j),用trapz函数求其和s2=trapz(s1);(3)s2(k)又是随x而变的,再沿x方向用trapz函数积分.由于MATLAB的数组运算可以代替一个for循环,所以三重积分只用了两组for语句.使本题的程序比较简明。 三重积分程序exn546dx=0.01;dy=dx;dz=dx;x=0:dx:1;fork=1:length(x)x1=(k-1)*dx;y=0:dy:x1;forj=1:length(y)y1=(j-1)*dy;z1=0:dz:x1*y1;%z1数组f=x1.*y1.^2.*z1.^3;%f(z1)s1(j)=trapz(f)*dz;%沿z1积分ends2(k)=trapz(s1)*dy;%沿y1积分end,s=trapz(s2)*dx%沿x1积分 六．微分方程的数值积分MATLAB中用来进行常微分方程数值积分的函数有好多种，例如ode23,ode45,…等，ode是常微分方程(ordinarydifferentialequation)的缩写。它们都用来解形如的一阶微分方程组在给定初始值y0时的解。对入门者而言，会一种最简单的ode23就行。它的最简单的调用格式为：[T,Y]=ODE23(ODEFUN,TSPAN,Y0)其中，输入变元TSPAN=[t0,tf]是自变量的初值和终值数组，Y0是输出变量向量的初值，ODEFUN则是描述导数的函数f(y,t)。很大一类微分方程都可以用这种一阶微分方程组（或向量形式的微分方程）描述，关键就是会列写导数函数f(y,t)，下面举例说明。 微分方程数值积分【例5-4-7】用数值积分法求解微分方程设初始时间t0=0;终止时间tf=3π;初始条件y(0)=1,y’(0)=0.解：先将方程化为两个一阶微分方程的方程组，其左端为两维变量的一阶导数。 微分方程化为标准形式写成矩阵形式为其中为取代变量y的变量向量，为x的导数，在程序中用xdot表示。x的初始条件为这就是待积分的微分方程组的标准形式。用MATLAB语句表述为：xdot=[0,1;-t,0]*x+[0;1]*(1-t^2/pi^2); 【例5-4-7】数值解的程序将微分方程的右端写成一个exn547f.m函数程序,内容如下：functionxdot=exn547f(t,x)u=1-(t.^2)/(pi^2);xdot=[0,1;-t,0]*x+[0;1]*u;%向量导数方程主程序exn547如下,它调用MATLAB中的现成的数值积分函数ode23进行积分。clf,t0=0;tf=3*pi;x0=[1;0];%给出初始值[t,x]=ode23("exn547f",[t0,tf],x0)%此处显示结果y=x(:,1);%y为x的第一列plot(t,y),grid%绘曲线xlabel("t"),ylabel("y(t)") 数值解程序exn547的运行结果◆程序运行的结果见图5-37。这个数值积分函数是按精度要求自动选择步长的。它的默认精度为1.e-3，因此图中的积分结果是可靠的。若要改变精度要求，可在调用命令中增加备选变元，具体做法可键入helpode23查找。 exn547的运行结果的讨论从物理意义看，这个方程表示了一个变系数的无阻尼振动方程。如果这是一个机械振动，则弹簧刚度随时间成正比地增强，振动频率随之逐渐提高。为了看得更为清楚，设弹簧刚度随时间按三次方增强，即方程的第二项系数为y3，则只要把子程序exn537f中的核心语句改为xdot=[0,1;-t.^3,0]*x+[0;1]*u;重新运行程序exn537，就得到频率迅速提高的波形,如图5-37。如果我们在原来的方程中加进y的一阶导数项（阻尼项），也只要在函数子程序中把矩阵的系数作些改动，马上就可以得出新的结果。由此可见，用计算机解题的极大优越性。 七．常系数线性微分方程的数值解第四章4.3.5节介绍了常系数线性微分方程用MATLAB求解的问题。其实这类方程是有解析解的，这个解析解取决于微分方程系数多项式的根。而四次以上多项式的求根却没有解析解，这就要依靠MATLAB用数值方法解决代数问题，这个函数称为residue，根据微分方程左端系数多项式a和右端系数多项式b，就可求根p和求留数r。[r,p]=residue(b,a)读者可以参阅4.3.5节的算例，并可参阅第七章中的机械振动和第九章中求系统响应的例题 八.符号数学求不定积分不定积分问题要用符号数学的公式推理功能来解决问题。而根据本书的指导思想和风格，主要强调数值计算和计算的道理，不把公式推理放在主要的地位。但是工作中如果遇到这种需要，还是应该利用符号数学的功能来解决，就算是查积分表也是应该会查的。所以也大致地介绍一下例题的类型和解法。 符号数学解不定积分例548(a)例5-4-8(a)。求不定积分解：因为是不定积分，不能用数值方法计算，只能用符号数学工具箱。程序为：>>symsx,y=x^2*atan(x),>>Z=int(y)得到y=x^2*atan(x)Z=1/3*x^3*atan(x)-1/6*x^2+1/6*log(x^2+1) 符号数学求不定积分例548(b)例5-4-8(b)。解下列积分，画出解的曲线。解：程序为>>symsx>>Y=int("(cos(x))^2+sin(x)")系统运行的结果为：Y=1/2*cos(x)*sin(x)+1/2*x-cos(x)+C代入边界条件，得：>>C=1-(1/2*cos(pi)*sin(pi)+1/2*pi-cos(pi))结果是C=-pi/2 符号数学求不定积分例548(c)例5-4-8(c)。求下列积分。解：程序为：>>symsx,Y=int(exp(-x^2))%不定积分：>>Y1=int(exp(-x^2),0,1)%定积分：运行结果为：Y=1/2*pi^(1/2)*erf(x)Y1=1/2*erf(1)*pi^(1/2) 【应用篇中与本节相关的例题】①【例6-5-1】和例6-5-2磁场计算由元电流生成的磁场积分，得到全部线圈产生的磁场。②【例7-1-3】导弹跟踪目标的轨迹，根据x,y两个方向的速度积分得出轨迹。③【例7-1-5】有空气阻力下的抛射体轨迹，是微分方程的数字积分求解问题。④【例7-2-2】懸臂梁的桡度计算，这是纯粹的两次积分问题。⑤【例7-2-3】简支梁的桡度计算，也是纯粹的两次积分问题。⑥【例7-3-1】，【例7-3-2】，一自由度自由振动和强迫振动，都是常微分方程求解的典型问题。这两道题用的是解析解，只是用数学软件来求根和绘制波形。 【应用篇中与本节相关的例题】⑦【例7-3-3】两自由度振动，是四阶的矩阵微分方程的问题。本题对于高阶矩阵指数采用了数值解，并且用矩阵对角化的方法，对振动的模态作了分解和分析。⑧【例9-1-2】任意高阶连续线性系统冲击响应的计算，它实际上是求常系数线性微分方程的解析解，可归结为代数问题，用MATLAB多项式函数库来求解。⑨【例9-1-4】多个放大器串联的脉冲响应仍然属于常系数线性微分方程在初始条件下的解，本题的特殊性在于遇到了重根。要在理论上解决这个问题，相当麻烦，但是用MATLAB作一下数值处理，就可以非常方便地求解。从这里可以看出工程师和数学家处理数学问题的区别。⑩【例9-4-1】在电磁场课程中常常遇到偏微分方程，本例给出了它的数值算法。这种方法对于电磁场的计算有一定的普遍意义。 5.5节线性代数 线性代数的主要特点代数问题要靠大量四则运算来解题，当用笔算或计算器算时，用的模型是单个数与数的运算。方程组的元数N愈高，运算的次数按N的平方或立方增加，出错的概率也迅速增长，所以笔算只能解三阶以下的问题。用这种方法去解高阶的问题，读者就感到抽象、冗繁和枯燥。用计算机算时，利用的是矩阵模型，它的运算对象是庞大的数据群组成的矩阵。只要给出原始数据列成的矩阵，键入一、两个命令就可得出准确的结果。所以它把冗繁变成了简明；MATLAB的作图能力很强，容易把抽象的问题形象化；又由于解题简捷，很容易在课程中引入并解决大量的应用例题，使得原本枯燥的课程变得丰富精彩。线性代数解决的实际问题大体上就是三类。①求线性代数方程组（包括欠定、适定和超定）的解；②分析向量的线性相关性；③矩阵的特征值与对角化。下面的例题主要围绕这几个方面展开。 消元法的MATLAB语句【例5-5-1】用矩阵的初等变换将消元为上三角阵，求等价变换乘子B,检验其行列式、秩和迹。解：程序exn551如下：A=[107;415;2-19];A0=A%输入A，留一备份A(2,:)=-A(2,1)/A(1,1)*A(1,:)+A(2,:);%消去A(2,1),A1=A,B1=A1/A0%求本次等价乘子B1A(3,:)=-A(3,1)/A(1,1)*A(1,:)+A(3,:);%消去A(3,1)A2=A,B2=A2/A1%求本次等价乘子B2A(3,:)=-A(3,2)/A(2,2)*A(2,:)+A(3,:);%消去A(3,2)A3=A,B3=A3/A2%求本次等价乘子B3B0=A3/A0%求总的等价乘子B0 从例5.5.1中要学到什么？请读者从三次消元中归纳出消元法的语法规则。如果选第i行为基准行，其第k列的元素为基准元素，要把第j行的第k列的元素消元为零，则应该执行如下运算：A(j,:)=-A(j,k)/A(i,k)*A(i,:)+A(j,:)可以专门编成一个消元子程序。读者还可观察这几个初等变换矩阵的构成特点。不难验证B0=B3*B2*B1。要注意，这几个乘子相乘的次序是不能颠倒的。 【例5-5-2】解下列代数方程组：解：可以把线性方程组写成矩阵方程A*x=b，其中:解这个矩阵方程可以用下列几种方法： 方法一最简行阶梯法用将其增广矩阵[A,b]化为最简行阶梯型。相应MATLAB命令为rref。程序如下：A=[6,1,6,-6;1,-1,9,9;-2,4,0,-4;4,2,7,-5];b=[7;5;-7;-9]U=rref([A,b])程序运行的结果为：这个矩阵代表了以下的等价方程组，所以等式右端的四个数也就是方程的解。 其他方法及比较方法二：用x=inv(A)*b方法三：用x=Ab在方程数与未知数数目相等的适定方程的情况下，三种方法所得的结果是一样的。如果方程是欠定的，则行列式det(A)≈0。方法二、三会得不出可信的解，如果方程是超定的，用方法二将导致出错警告，用方法三将得出最小二乘意义下的解。为了在任何条件下得出可靠的结果，建议经常用化为最简行阶梯型的方法来解方程。 欠定方程解例【例5-5-3】把例5-5-2中第四个方程改为：，求其解。解：输入新参数A=[6,1,6,-6;1,-1,9,9;-2,4,0,-4;4,2,7,-778/222];b=[7;5;-7;877/222];方法一：键入U=rref([A,b])，得到这个最简行阶梯形式说明原来的方程组是欠定的。 欠定方程组解的特点它等价于下列方程组：这是一组包括四个变量的三个有效方程。因此没有唯一的解。其中x4可以任意设定，即可以看作任意常数c，：代入不同的c可以得到不同的解，因此欠定方程组有无数个解。这些解组成一根在空间中的直线。 用方法二、三：键入x=inv(A)*b或x=Ab，得到：Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=3.590822e-018.及计算机告诉我们，这个结果是不准确的。其原因在于矩阵A的行列式接近于零，不难检验，det(A)=0，rank(A)=3。就是说，A的逆极小。从逆条件数RCOND=3.590822e-018，可以得知计算结果有效数位减小的程度。MATLAB的是有效数位是16位，现在算得的结果有效数位要减去18位，所以得出的结果是根本不能相信的。 不相容方程组的示例如果把第四个方程的右端常数仍取为-9，则其行阶梯变换的结果为：最后一个方程成了一个矛盾方程0=1。这说明方程组不相容，无解。由此也可以看出，线性方程组求解最好还是用行阶梯简化的方法。因为它可以给出线性方程组的特征，避免计算的盲目性。 【例5-5-4】多项式插值问题给出平面上4个点的坐标值如右表，(1)。求对它进行插值的三次多项式，(2)。求t=1.5处f的近似值。(3)。如果要求此多项式多通过一点(-1,5)，求其系数。解：用多项式来插值，令它在四点上的值与表中相同，，得到ti0123f(ti)30-16 多项式插值问题的解(1)这个矩阵方程达到四阶，应该用计算机辅助求解了，编出MATLAB程序exn554C=[1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27],b=[3;0;-1;6],U=rref([C,b])得到，多项式为 多项式插值问题的解(2)（2）把t1=1.5代入此多项式，键入：f1=3-2*1.5-2*1.5^2+1.5^3得到：f1=-1.125。可以用以下语句画出插值图形，ezplot("t^3-2*t^2-2*t+3",[-1,4]),gridon,holdonplot(0:3,[3,0,-1,6],"*")plot(1.5,-1.125,"or")可以得到图5-39，曲线通过了图中四个给定的插值点（用*号表示），圆圈为f(1.5)的位置。 多项式插值问题的解(3)(3)。若要此三次多项式多通过一点(-1,5)，将此点坐标代入后得，方程组就成为五个方程四个未知数，很可能是矛盾的，要靠计算来判断。用以下程序：A1=[1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27;1,-1,1,-1],b1=[3;0;-1;6;5],U01=rref([A1,b1])得到 多项式插值问题的解(4)注意到系数增广矩阵最后一列是常数项，得出右方的同解方程组。其最后一个方程成为0=1,说明方程组是不相容的，属于超定方程，无通常意义下的解。用MATLAB的pinv函数可以求出它的最小二乘解，键入：a1=pinv(A1)*b1，得：多项式成为 多项式插值问题的解(5)画出它的图形，可以继续键入t=-1:0.1:4;plot(t,a1(1)+a1(2)*t+a1(3)*t.^2+a1(4)*t.^3,":r")plot(-1,5,"x")此曲线也用虚线画在图上。它不通过给定的5个点中的任何一个，说明都有误差。误差向量E及五项误差的平方和EE可以用矩阵方程计算如下E=A1*a1-b1，EE=norm(E)求得，结果为 【例5-5-5】化学方程的配平配平下列小苏打与柠檬酸的反应。解：建模，按每种物质的元素Na,H,C,O写成列向量，按四种元素左右平衡列出四个方程，得： 化学方程配平程序exn555：>>A=[1,0,-3,0,0;1,8,0,-2,0;1,6,-6,0,-1;3,8,-7,-1,-2],b=[0;0;0;0]>>U=rref([A,b])结果为整数格式 配平程序exn555运行结果由此可见，x5必须最小取30，才能保证各个分量xk为整数，于是得到最后可以得出配平后的化学反应方程为这个待配平的化学方程有四个方程和五个未知数，所以它是一个欠定的线性方程组，其解乘以任何常数仍然为方程的解。不过当我们另加了一个条件——系数取最小的正整数，结果就是唯一的了。 【例5-5-6】向量的线性相关性【例5-5-6】设由三个在R4空间的四维基向量v1,v2,v3张成的子空间，问w1和w2是否在此子空间内。其中：解：本题的要点在于研究w是否能由v1,v2,v3能以线性组合的方式组成，即是否找到三个常数c1,c2,c3，以便得到三个基向量的线性组合：能实现这个关系的w，就称为与向量组v1,v2,v3线性相关。反之就线性无关。 方法一用最简行阶梯型◆从向量组的最大无关组的概念，可以解决这个问题。把五个向量列成一个矩阵A，对它作行阶梯变换。键入：A=[v1,v2,v3,w1,w2],U0=rref(A)得到：可见最大无关组是由v1,v2,v3,w1四个向量组成的，即这它们线性无关，w1不可能由v1,v2,v3组合而成。w2与v1,v2,v3构成的矩阵却只有三行，所以它与向量组v1,v2,v3线性相关。U0矩阵中对应于w2列的系数恰好就是c。即c1=-1,c2=-2,c3=1。 方法二用向量组的秩判别把三个列向量并排成矩阵v[v1,v2,v3]，c1,c2,c3排列成列向量c，则上述方程可以写成矩阵与向量相乘的形式v*c=w。若w与v线性相关，其组合矩阵[v,w]的秩应该与v的秩相同，反之，其秩应该加1。故列出程序ag556：v1=[7;-4;-2;9];v2=[-4;5;-1;-7];%输入参数v3=[9;4;4;-7];w1=[-9;7;1;-4];w2=[10;-2;8;-2];v=[v1,v2,v3];%将三个基向量组成矩阵dr1=rank([v,w1])-rank(v)%v与w1组合后矩阵秩的增量dr2=rank([v,w2])-rank(v)%v与w2组合后矩阵秩的增量运行这个程序，得到dr1=1，dr2=0这说明w1不是v1,v2,v3的线性组合，而w2则是v1,v2,v3的线性组合，w2将位于v1,v2,v3所张成的R3子空间内。 矩阵的条件数和计算精度(1)【例5-5-7】三阶hilbert矩阵H3=hilb(3)如下规律构成：现要构成一个四元一阶线性方程组H4*x=b4,要求（1）求H4的行列式D=det(H3)及条件数c=cond(H4)；（2）设x的四个分量均为1，求b4；（3）按方程H4*x1=b4求x1,并与x相比较，求出其最大相对误差。 矩阵的条件数和计算精度(2)解：程序exn557如下：H4=hilb(4),D=det(H4),c=cond(H4)x=ones(4,1);b4=H4*x,formatlong,x1=H4b4,e=max(abs(x1-x)./x)程序运行结果为：H4=1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000.16670.25000.20000.16670.1429D=1.6534e-007c=1.5514e+004e=2.958744360626042e-013 矩阵的条件数和计算精度(3)H4的行列式D约为10-7，相当接近于零了。方程是否该算成是奇异的？这不是由数学家所能定的，而是取决于工程的判据。这时主要看条件数c。c达10的四次幂以上，意味着解的有效数位将减少四位。由于MATLAB的计算精度达到16位有效数位，故解的相对误差仍可保持在10–12以下（算出的e可以证明），这在工程上是完全可以接受的。读者可自行扩展此程序，把H矩阵的阶次逐步提高，看到什么阶次相对误差会达到百分之百。到什么阶次MATLAB会发出警告。 求方阵的特征值和特征向量【例5-5-8】设有对称实矩阵，试求其特征方程、特征根和特征向量，并讨论矩阵的行列式和迹与特征值的关系。分三个步骤：不用eig函数，按照原理，利用MATLAB函数分步求解；用eig函数求解，进行对照校核；求特征值矩阵的行列式与迹。 例5-5-8的MATLAB程序exn558%（1）分步计算特征值和特征向量A=[2,4,9;4,2,4;9,4,18],%输入矩阵参数f=poly(A)%求其特征多项式的系数向量fr=roots(f)%求其特征根fori=1:3%将三个特征根分别代入特征方程p(:,i)=null(A-r(i)*eye(3));%求齐次方程的基础解end,p,d=diag(r)%列出特征向量矩阵和特征根矩阵%（2）用eig函数求特征值和特征向量：[p1,d1]=eig(A)%一步求出特征值和特征向量%求原矩阵及特征根矩阵的行列式和迹det(A),det(d),trace(A),trace(d), (1)分步计算的结果(1)分步计算的结果为：f=1.0000-22.0000-37.0000122.0000r=23.3603-3.06451.7042p=0.4149-0.84830.32900.24190.45140.85890.87710.2767-0.3925说明此矩阵A的特征方程为，特征根为23.3603，-3.0645，1.7042。将特征根在主对角线上排列即构成特征值矩阵。 (2,3)用eig函数计算的结果为：(2)p1=-0.8483-0.32900.41490.4514-0.85890.24190.27670.39250.8771d1=-3.06450001.704200023.3603两者的差别仅仅是特征值和特征向量的排列不同，因为eig函数中把特征值按递增排列。（3）det(A)=-122,det(d)=-122.0000,trace(A)=22,trace(d)=22,可见其特征根的和仍为原矩阵的迹，其特征根的积仍为原矩阵的行列式。 【例5-5-9】对称实矩阵对角化设有对称实矩阵试求解：建模：矩阵指数只有在A是对角矩阵时才可按对角元素直接计算，即：注意它的非对角元素不符合指数运算规则，。 实矩阵对角化求矩阵指数(1)要利用这个性质，首先把A对角化。使A=pDp-1，将此式代入y的表示式中，可得：为求其特征根和特征向量矩阵D,p。编成程序exn556：A=[-2,4,1;4,-2,-1;1,-1,-3],[p,D]=eig(A)得p=-0.6572-0.26100.70710.65720.26100.70710.3690-0.92940.0000D=-6.5616000-2.43840002.0000 实矩阵对角化求矩阵指数(2)由此得到：算式中的最后一步可以用下列语句来实现。y=p*[exp(D(1,1)),0,0;0,exp(D(2,2)),0;0,0,exp(D(3,3))]*inv(p)得到 应用篇中与本节相关的实例(1)①【例6-1-1】的物理数据的处理，用一根直线去拟合有测量误差的数据点，这就是解一个超定的线性方程组。每个方程都有误差，找到的最佳结果是使各方程的误差均方值为最小，也就是最小二乘问题。②【例6-2-3】和【例7-1-2】的静力学平衡问题，因为涉及两个物体，平衡方程和待求变量会超过四个。这些通常都是线性方程组，用矩阵来解可以一次解出全部未知数，最为方便快捷。③【例7-2-1】的材料力学静不定问题和上述静力学平衡问题相仿，只不过加了几个变形协调方程，这些方程仍然是线性方程，所以不管方程和变量增加了多少，都只要用一个矩阵方程一次就解出全部结果。只要系数输入没有错，结果肯定是对的。 应用篇中与本节相关的实例(2)④【例7-3-3】的二多自由度机械振动问题更是用矩阵解高阶微分方程的典型。它不仅涉及代数方程的解，而且涉及特征值和特征向量的计算问题。⑤【例8-1-1】的直流电路分析是典型的线性代数方程组，【例8-1-5】的交流电路分析也化成线性代数方程组，只不过方程的系数矩阵和未知数向量都是复数。在解复数方程组时，MATLAB更显示了比手工解题的极大优越性，包括绘制相量图。【例8-1-8】网络参数都是用矩阵表示和矩阵运算的。⑥【例8-2-2】的晶体管放大电路分析得出了一个复数矩阵方程，它的所有系数元素又都是频率的函数。给了不同的频率就可以得出在该频率上的结果，设50个频点，用一个循环语句就可以快速解50个复数矩阵方程，算出50个点上的输出。这比手工计算提高效率将达千百倍之多。⑦【例9-1-3】高阶微分方程的零输入响应，这在数学上高阶线性微分方程的初值问题。求它的各个输出分量归结为解一个同阶的矩阵方程，其系数矩阵就是范得蒙特矩阵，在这个例题中作了公式推导并给出了数值解例。 应用篇中与本节相关的实例(3)⑧【例9-2-2】给出了离散傅立叶变换的程序，由时域采样信号求出它的频谱，可以看作一个复数矩阵与时间样本的乘积，它就是著名的傅立叶变换。这个变换矩阵是N×N的方阵。在实际应用中，N通常至少为1024。所以，计算量极大，MATLAB是这个行业必不可少的工具。⑨【例9-3-1】和【9-3-2】给出了用线性代数推导线性系统传递函数的一般方法，它大大优于传统教材上介绍的梅森公式。它的优点之一是理论上有严格证明，而梅森公式通常都不证的；它的优点之二是完全可以靠计算机自动完成，所以复杂系统的传递函数推导必须依靠矩阵和计算机结合的方法。以上所举出的例子是不完全的。由于只是最近十多年微计算机才成为科技人员的手头工具，矩阵方程刚刚才变得简便易解，各门后续课程的教材许多方面还习惯地用老办法解题，该用而没有用矩阵，应该还有更多的实例。 5.6节概率论与数理统计 一各种统计分布函数表5-6-1中列出了20种分布函数，统计工具箱提供了包括分布函数~cdf(CumulattiveDistributionFunction)、概率密度函数~pdf(ProbabilityDistributionFunction)、分布函数的逆函数~inv(InverseCumulattiveDistributionFunction)以及这些分布的理论统计特性（均值和方差）计算函数~stat。 表5-6-1MATLAB中表示各种概率分布的前缀文字连续（数据）连续（统计量）离散（数据）贝塔分布(beta~)χ2分布(Chi2~)二项式分布(bino~).指数分布(exp~)非中心χ2分布(ncx2~)离散均匀分布(unid~)Г-分布(gam~)F-分布(f~).几何分布(geo~).对数正态分布(logn~)非中心F-分布(ncf~)超几何分布(hyge~)正态分布(norm~)T-分布(t~).负二项式分布(nbin~)瑞利分布(rayl~)非中心T分布-(nct~)泊松分布(poiss~)均匀分布(unif~)韦伯.分布(weib~) 表5-6-1的使用方法把表中不同分布后面括号中的文字作为前缀，把所需计算的特性作为后缀，就可以组合成一个函数。例如离散类的二项式分布有binopdf,binocdf,binoinv,binostat四个函数，连续类数据的标准正态分布有normpdf,normcdf,norminv,normstat四个函数，连续类统计量的χ2分布有chi2pdf,chi2cdf,chi2inv,chi2stat四个函数，…等等，20种分布就有80个函数。由于各种分布函数的解析形式都是已知的，这些子程序的编写并不困难。 【例5-6-1】概率分布曲线的绘制【例5-6-1】(a)求标准正态分布N(0,1)，自由度V=10的χ2分布和N=10,p=0.2的二项式分布B(10,0.2)的分布函数，并画出其概率密度曲线和分布曲线；(b)。求出分布函数为0.05和0.95处的随机变量值；(c)。求出这几个分布的均值和方差。解：分别键入helpnormpdf,helpbinopdf,helpchi2pdf以了解它们的用法，特别是了解输入变元的意义，得知：f=normpdf(x,mu,sigma)其中x为随机变量数组，mu为均值，sigma为标准差。f=chi2pdf(x,V)其中x为随机变量数组，V为自由度数（整数），f=binopdf(x,N,p)其中x为随机变量整数数组，N为次数，0≤p≤1为成功概率 概率分布参数的选择题目中给出的变元参数应足以用来调用这些概率密度函数，只有随机变量x的取值及范围，需要事先对该分布的特性有所了解。首先要弄清它是离散量还是连续量，其次要取适当的范围。范围取小了不能显示分布的全局，取大了又可能显示不出细节。正确的取法应该使该范围内分布函数F(x)的左边界值略小于0.025，右边界值略大于0.975，这就可以基本涵盖随机变量以概率95%存在的主要区域，又不致涉及关系不大的区域。初学者往往要作几次试探才能做到。好在在计算机上改几个参数、作几次试探是很简便的事。 绘制概率分布函数的程序exn561x1=[-3:0.1:3];%标准正态随机变量取值范围f1=normpdf(x1,0,1);%标准正态分布的概率密度函数F1=normcdf(x1,0,1);%标准正态分布的分布函数subplot(2,2,1),plot(x1,f1,":",x1,F1,"-"),gridon%绘曲线line([-4,4],[0.025,0.025]),line([-4,4],[0.975,0.975])%画横线x2=[0:0.1:20];%试探得出的范围f2=chi2pdf(x2,10);%χ2分布的概率密度函数F2=chi2cdf(x2,10);%χ2分布的分布函数subplot(2,2,2),plot(x2,f2,":",x2,F2,"-"),gridon%绘曲线line([0,20],[0.025,0.025]),line([0,20],[0.975,0.975])%画横线);…….运行程序得出图5-40，其中实线是分布函数的曲线，虚线则是密度函数的曲线。上下两根横线分别为0.975和0.025， 三种分布的概率密度和分布曲线 对分布函数图形的改进第三个子图中，随机变量是离散取值的，在两个相邻的取值点之间的概率不会变化，所以它的分布函数表现为阶梯形。密度函数是分布函数的导数，所以在阶梯突变处导数为无穷大，宽度又为无穷小，面积等于阶梯的高度，通常用一个脉冲表示。脉冲的高度表示它所包含的面积，也就等于阶梯波形的高度。而plot函数画图是把各离散点之间用直线联接，所得的曲线是不对的。应该要用stairs命令画阶梯图，用stem命令画脉冲图。改正后的程序如下：subplot(2,2,4),stairs(x3,F3,"-"),stem(x3,f3,":"),图中第四子图给出了改正后的结果，可见其概率密度函数是离散的脉冲。从中还可以判断，x3的取值范围不必为[0:10]，取[0:5]就够了。另外，第二子图上的密度函数波形太小，如果对f和F取不同的纵坐标，那样可以得出更好的图形。 分布函数逆函数的用途(b)给定分布函数F=α（α1）求出的x，简称α下分位点，习惯上用λα表示。这和已知x求cdf恰好是逆函数的关系，即输入变元与输出变元恰好调换了位置，对正态分布情况，逆函数norminv的调用格式为X=NORMINV(F,MU,SIGMA)其中F为给定的分布函数值，而X为对应的随机变量边界值。题目中给定了的两个边界值Fb=[0.05,0.95]，即求对应的随机变量x的边界值xb。，随机变量在上下两个边界值xb(即[λα/2，λ1-α/2])之间取值的置信度等于1-α，其他参数与normpdf和normcdf中的相同。对于其他分布，可依此类推，不再赘述。 用逆函数求上下分位点将边界值Fb作为输入变元，可求得相应的分位点。：alpha=0.1;%取双边90%的置信区间Fb=[alpha/2,1-alpha/2]%输入变元%三种分布的双侧α分位点lambda1=norminv(Fb,0,1)lambda2=chi2inv(Fb,10)lambda3=binoinv(Fb,10,0.2)程序运行后得到lambda1=-1.64491.6449lambda2=3.940318.3070lambda3=04这就是三种分布函数下的双边90%置信区间改变alpha的值，可得到各种不同置信度下的置信区间。 各种分布理论统计特性的计算(c)。对以上3种分布，MATLAB还给出了它们的理论统计参数，即理论均值Mu和方差值var的计算方法，所用的命令为~stat。例如：[Mu1,var1]=normstat(0,1)[Mu2,var2]=chi2stat(10)[Mu3,var3]=binostat(10,.2)运行结果为：Mu1=0，var1=1Mu2=10，var2=20Mu3=2，var3=1.6000本题虽然只给出了3种分布的计算，这些概念和方法可以类推到表5-6-1中20种分布的理论计算中。 概率分布演示工具disttool键入disttool就会出现一个图形界面(右图)，其中有分布曲线，周围有各种参数和类型的选择窗。可以很方便地用鼠标操作改变参数和选择类型，得到相应的曲线。 二。随机样本数据的生成stats工具箱同时也提供了这20种分布的随机数生成程序~rnd。均匀分布随机数的计算机生成本身就是一个研究了几十年的专题，其他分布的随机数通常又由均匀分布随机数进行变换而得。本书将不讨论它们的编程，只着重于它们的应用。下面的例子将说明这类函数的调用方法。【例5-6-2】按例5-6-1的三种分布，分别生成10000个随机数，画出它们的直方图（分布图），并计算各自的数学期望和方差。 随机数生成及其直方图绘制各种分布的随机数生成函数名为~rnd。其中表示分布类型的前缀~由表5-6-1给出。绘制统计数据Y直方图的命令为hist(Y,N)，其中N为直方图的分区数目，其缺省值为10。例如：Y1=normrnd(0,1,1,10000);subplot(2,2,1),hist(Y1,50)Y2=chi2rnd(10,1,10000);subplot(2,2,2),hist(Y2,50)Y3=binornd(10,0.2,1,10000);subplot(2,2,3),hist(Y3,50)得出的图5-42。 三种分布的随机数统计分布图 实际统计参数的计算其中第三个分图的直方图各条宽度不均匀，这是由于二项式分布是离散数据，它的取值最大只可能到10。在目前实际数据中，最大只到8。因此，可改用以下的命令：subplot(2,2,4),hist(Y3,8)得到符合一般直方图要求的第四个分图。把这些图与图5-40中的密度分布曲线相比，可见是非常相似的。实际随机变量的样本均值Xbar,样本的标准差sigma和样本方差s2(标准差的平方)可用以下命令计算Xbar1=mean(Y1),sigma1=std(Y1),s21=var(Y1)得出Xbar1=0.0246,sigma1=1.0059,s21=1.0119Xbar2=mean(Y2),sigma2=std(Y2),s22=var(Y2)得出Xbar2=9.9947,sigma2=4.4453,s22=19.7605Xbar3=mean(Y3),sigma3=std(Y3),s23=var(Y3)得出Xbar3=1.9745,sigma3=1.2520,s23=1.5676 演示工具randtool生成的图形界面键入randtool也会出现一个图形界面，如右图。其中有实际生成的随机数样本分布的曲线，周围有各种参数和类型的选择窗。可以方便地用鼠标和键盘改变参数，得到相应的曲线。 三。用样本推断总体的统计量由实际样本值算出的统计量与例5-6-1中的理论值之间存在误差，而且样本统计量也是随机量。计算机每运行一次，所得到一组观察值Yi的统计量也会不同。它们也有各自的分布规律，比如对于μ,σ2相同的正态分布总体，各样本组的标准化均值按Z-分布，各样本组的标准化方差按χ2分布。当样本数取得很大时，根据中心极限定理，分布的误差就很小。这样得出的是总体统计量的一个（有误差的）估值，称为点估计法。如果样本数小，点估计法误差就会很大，甚至对人们造成误导。这时较好的方法是给出总体统计量一个取值区间（置信区间），并给出它在此区间取值的概率（置信度）。在实际工程中，做试验往往是很费时费钱的。通常希望用最小的试验样本集来推断总体的统计量所取的数值范围，并且希望有相当的精度和可信度。 正态分布样本与总体统计量的关系(a)。已知总体方差σ2，但均值μ未知。要由样本均值和方差S2推断总体均值μ。这时，样本均值标准化后的Z值满足正态分布：（5.6.1）(b)。总体均值μ和方差σ2均未知，要由样本方差S2推断总体方差σ2。这时，样本方差S2标准化值满足χ2分布：（5.6.2）(c)。总体均值μ和方差σ2均未知，要由样本均值和样本方差S2估计总体均值μ。μ经过样本均值和样本方差S2标准化后满足T-分布：（5.6.3） 【例5-6-3】对某种产品的尺寸进行检验，由长期的统计知道其总体方差为σ2=4（微米2），当日抽测其中10个工件的偏差分别为2;1;-2;3;2;4;-2;5;3;4（微米），试对零件长度偏差的平均值（总体的数学期望）作出置信度为90％的区间估计。解：在MATLAB程序中，用Xbar表示均值，s2表示S2，则样本均值和样本方差分别为：x=[2;1;-2;3;2;4;-2;5;3;4];Xbar=mean(x),s2=var(x)得到：Xbar=2,s2=5.7778因为对于正态分布的随机样本，其样本均值Xbar和总体均值μ误差的标准化结果Z也服从正态分布。要估计总体均值μ的范围。要先找到其90％置信区间边界点[λ0.05，λ0.95]，这可用正态分布逆函数查出，其语句为：lambdam=norminv([0.05,1-0.05]) 【例5-6-3a】总体均值μ的区间估值根据Z的这样两个边界，可以计算相应的μ的估值区间。把（5.6.1）式进行移项，得到，把公式中的Z用它的两个置信区间边界点来置换，在MATLAB中令varx=σ2，并利用数组计算方法，同时计算两个边界的值，即把公式表达为：相应的MATLAB程序exn563就是：x=[2;1;-2;3;2;4;-2;5;3;4];Xbar=mean(x),s2=var(x)lambdam=norminv([0.05,0.95])n=length(x),varx=4Mut=Xbar-lambdam*sqrt(varx/n)%Mut表示总体均值的估计范围 程序的运行结果为：lambdam=-1.64491.6449n=10varx=4Mut=3.04030.9597于是根据这十个样本，我们有90%的把握来断定，该工件尺寸的总体平均值将大于标称值0.9597～3.0403（微米）。如果我们只取一个误差为+2的样本，其他参数都不变，这时算出的置信区间将为：Mut=5.2898-1.2898这意味着，我们甚至无法估计工件的尺寸总体值是正偏差还是负偏差。所以取的样本太少，是很难有信心较准确地推断工件的总体偏差的。 例5-6-4最大飞行速度方差的估计飞机的最大飞行速度X~N(μ,σ2)，但μ,σ2未知。现对飞机的速度进行14次独立测试，测得数据如下(单位：m/s)：422,419,426,420,426,423,432,428,438,434,412,417,414,441;试按置信度0.95对总体X的均值μ和方差σ2进行区间估计。解：先求标准差的区间估值，总体均值μ和方差σ2均未知，要由样本方差S2推断总体方差σ2，属于情况(b)，方差按χ2分布。其区间估值可以分成两个步骤：(1)由，找到的置信区间边界点，(2)把上式移项整理为：求得σ2的估值区间。 求总体方差估值的程序exn564先用MATLAB语句输入数据，然后进行这两个步骤。在程序中设样本方差为s2，样本标准差为sigma，总体标准差为sigmat。所对应的MATLAB程序exn564如下x=[422,419,426,420,426,423,432,428,438,434,412,417,414,441]";n=length(x),s2=var(x),Xbar=mean(x)%根据给出的α,反查χ2的置信区间边界点。因已知置信度为0.95,故α应取1-0.95=0.05alpha=input("alpha="),%由α查χ2置信区间lambdach=chi2inv([alpha/2,1-alpha/2],n-1)%由置信区间求估值区间sigmat=sqrt((n-1)*s2./lambdach) 程序exn564运行结果程序运行时，若按其提问‘alpha=’，键入0.05，则所得的结果为：n=14s2=76.4396Xbar=425.1429alpha=0.0500lambdach=5.008824.7356sigmat=14.08536.3383再求总体均值的估值区间，因为我们对于总体方差并不知道，而要靠样本数据来估计总体均值μ。这就属于上述情况(c)。 例5-6-4最大飞行速度均值的估计解：由找到的置信区间边界点，再由求出μ相应的估值区间。这两个步骤所对应的语句为%由α查T-置信区间lambdat=tinv([alpha/2,1-alpha/2],n-1)%由置信区间求估值区间Mut=Xbar-sqrt(s2/n)*lambdat运行的结果为：Mut=430.1909420.0948可见最大速度按95％的置信度应在420～430【米/秒】之间。 用工具箱函数进行估计用MATLAB提供的normfit函数来根据样本集的数据x直接估计出总体统计量。其调用格式为：[Xbar,sigma,mut,sigmat]=normfit(x,alpha)在本题中，用测得的x为输入变元，键入：[Xbar,sigma,mut,sigmat]=normfit(x,0.05)即可得到样本的统计量为：Xbar=425.1429sigma=8.7898（注意，即sigma=sqrt(s2)）总体的估计量为mut=420.0678430.2180sigmat=6.372214.1608如取置信度为90%，则用下列语句[Xbar,sigma,mut,sigmat]=normfit(x,0.1) 四。蒙特-卡罗随机实验法【例5-6-5】用随机实验法求圆周率估计值。解：先生成在-1≤(x,y)≤1中均匀分布的随机数作为坐标系（x,y）中x，y的值，这样两个坐标就决定图5-44所示的正方形中一个点，如果随机数的分布是均匀的，则落入圆中和方形区域中点的数目之比，就代表了圆面积和正方形面积之比。又已知正方形面积为4，这样就可以求出圆的面积，这个面积就等于圆周率的估计值。 随机实验求pi程序exn565N=input("取的总点数=")；%生成均值为零，在-1~1间均匀分布的随机数x,y%点(x,y)全在方形区域内,k为落入圆内部的点数，初值为零x=2*(rand(1,N)-0.5);y=2*(rand(1,N)-0.5);k=0;fori=1:Nifx(i)*x(i)+y(i)*y(i)<=1k=k+1;end%落入圆的中点数endpihat=4*k/N%圆面积＝正方形面积×k/N运行此程序，分别在程序的提示下键入N=10000,100000,1000000，...等，可以得出圆周率的估计值pihat为3.1512，3.1448，3.1417,...，输入点数愈多，得出圆周率的精度就愈高. 五。统计工具箱中的其他函数前面介绍的只是MATLAB的统计工具箱(stats)中的一些初级的功能。此工具箱还有大量的功能很强的函数。这些功能包括聚类分析，回归分析，假设检验，多变量统计分析，试验设计，统计过程控制等等。这些内容的应用性很强，MATLAB针对各种应用编写了相应的函数，往往每一个应用有一个专门的函数,只讲输入输出变元的意义和用法,不讲原理，现举一个假设检验的例子说明。 假设检验函数ttest的用法(1)【例5.6.6】(a)用例5.6.3的数据，检验该样本集的总体均值是否可能为零？(b)用例5.6.4的数据，检验该样本集的总体均值是否可能为425？取置信度为0.9。解：因为总体方差为未知，故总体均值服从T-分布，MATLAB为这样的分布进行总体均值估值的假设检验编写的函数为ttest。键入helpttest可以得到它的用法：TTEST假设检验:将样本均值和一个常数进行T检验，调用格式为[H,P,CI,STATS]=TTEST(X,M,ALPHA,TAIL)判断满足正态分布的样本X是否可能具有总体均值M。其他输入变元的意义：ALPHA为1-置信度，TAIL是用来对H=1的意义进行控制的参数（见输出变元部分），默认值为M=0,ALPHA=0.05.，TAIL=0。 假设检验函数ttest的用法(2)输出变元的意义如下H=0表示：“总体均值=M”的假设成立；H=1表示：“总体均值=M”的假设不成立；其进一步的意义，则随输入变元TAIL的不同，意义为若TAIL＝0，表示：“总体均值M”；若TAIL＝1，表示：“总体均值>M”；若TAIL＝-1，表示：“总体均值