第十六讲在高等数学中的应用

第十六讲在高等数学中的应用数值解还是解析解？在“语言篇”介绍了MATLAB数学软件基本语法的基础上，本篇讨论如何用它来解决大学所学过的数学问题。本篇的第1～4节依次讨论微积分中的极限与导数、解析几何、数列和级数以及数值积分问题。主要通过一些例题说明如何灵活使用MATLAB的各种函数来解题。对于学生而言，各门后续课程和未来的工程实践中遇到最大量的将是数值计算问题。计算机首先是计算的工具。计算机的计算过程和方法都是从计算器升级而来的，学生可以理解和接受其每一步，甚至自己都可以编出相应的程序，这是数值计算的一个长处。其次，用推理方法只能解决很少一部分有解析解的数学命题。比如许多函数是无法求不定积分的，而它们的数值积分却都可以求得。因此优先让学生掌握数值方法等于给他们教会了具有更普遍适用的科学计算方法。对于他们今后的工程生涯将有更广泛的用处。当一个例题可以同时用数值方法和符号推理方法解决时，我们采取的原则是优先讲数值方法。因为通常这两种方法在编程上很相似，但又有一些关键性的差别，初学者很容易混淆。对于这类读者，还是先掌握数值方法为好。5.1函数、极限和导数一．单变量函数值的计算和绘图【例5-1-1】单变量函数的计算和绘图,设要求以0.01s为间隔，求出y的151个点，并求出其导数的值和曲线。解:43sin4323tyetdt=0.01;t=0:dt:1.5;%设定自变量数组tw=4*sqrt(3);%固定频率y=sqrt(3)/2*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)/dt;subplot(2,1,2)plot(t(1:length(t)-1),Dy)gridylabel('Dy(t)')%加标注set(gcf,'color','w')参变方程的计算和绘图【例5-1-3】摆线的绘制当圆轮在平面上滚动时，轮上任一点所画出的轨迹称为摆线。如果这一点不在圆周上而在圆内，则生成内摆线；如果该点在圆外，即离圆心距离大于半径，则生成外摆线。摆线绘制的程序◆建模：''sin()cos()(sin())cos()sin()cos()AAAAArrtirrtirtjrrrrtrtirtjxrtrtyrt00:1;1;0sin:cosAArxttyt令则◆MATLAB程序exn512：t=0:0.1:15;%设定参数数组x=t+sin(t);y=cos(t);%计算x,yplot(x,y)axis('equal')%绘图grid动画t=0:0.001:15;%设定参数数组x=t+sin(t);y=cos(t);%计算x,ycomet(x,y)摆线绘制程序的结果设r=1,令R=r,R=0.7及R=1.5时得到的摆线、内摆线和外摆线都绘于图5-1-3中。为了显示摆线的正确形状，x,y坐标保持等比例是很重要的，因此程序中要加axis(‘equal’)语句。02468-101234摆摆摆摆摆摆摆摆三曲线族的绘制【例5-1-4】三次曲线的方程为,试探讨参数a和c对其图形的影响.图5-1-3c和a取不同值时y=ax3+cx的曲线族解:◆方法因为函数比较简单,可以直接写入绘图语句中,用循环语句来改变参数.注意坐标的设定方法,以得到适于观察的图形。给出的程序不是唯一的,例如也可用fplot函数等,3yaxcx画曲线族的程序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-4程序运行的结果-202-3-2-10123xyc=0c=-3c=3-202-3-2-10123a=-3a=3a=05.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数列用for循环的表示方法计算n!（n的阶乘），它应该写成prod(1:n)fork=1:6x(k)=prod(1:k);endx三.函数项级数【例5-3-2】利用幂级数计算指数函数解：◆原理：指数可以展开为幂级数%输入原始数据，初始化yx=input('x=');n=input('n=');y=1;%将通项循环相加n次，得yfori=1:ny=y+x^i/prod(1:i);endy分别代入x=1,2,4,-4四个数，取n=10，结果如下表21......2!!nxxxexn用级数计算指数的结果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例用付利叶级数合成方波解：◆建模设方波的周期为2π,半周期为π.幅度为正负0.7854,如右图所示.它应该可以用傅立叶级数合成.用上题的程序分析,应有:55sin33sinsinttty现取不同的阶次,编写程序,画出其合成波形,在图上进行比较.MATLAB程序exn547◆设计如下：t=linspace(-pi,pi,201);%设定时间数组y=sin(t);plot(t,y),pausey=sin(t)+sin(3*t)/3;plot(t,y),pause%叠加三次谐波%用1，3，5，7，9次谐波叠加y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;plot(t,y)k=input(‘k=’)%K级谐波t=linspace(-pi,pi,201);%设定时间数组y=0forn=1:2:ky=y+sin(n*t)/n;endplot(t,y),5.4节数值积分和微分方程数值解求任意非线性方程f(x)=0的解【例】求下列方程的近似数值解。32102sin500yxxx解:可先大致看一下曲线的形状。它有三个过零点。分别在-10,-3,3附近。x=-12:0.01:5;y=x.^3+10*x.^2-2*sin(x)-50;plot(x,y)y1=zeros(1,length(x));holdonplot(x,y1,'k')grid-12-10-8-6-4-20246-400-300-200-1000100200300400实际上可直接用fzero求f(x)的根，方法是把f(x)定义为一个子程序函数ex517f.m，其内容为functiony=ex517f(x)y=x.^3+10*x.^2-2*sin(x)-50;然后执行x=fzero(‘ex517f’,x0)，即可得出靠x0较近的数值解x。二．数值定积分求面积【例5-3-1】用数值积分法求由x=0与x=10围成的图形面积，并讨论步长和积分方法对精度的影响。解:◆原理用矩形法和梯形法分别求数值积分并作比较，步长的变化用循环语句实现。MATLAB中的定积分有专门的函数QUAD,QUADL等实现。为了弄清原理，我们先用直接编程的方法来计算，然后再介绍定积分函数及其调用方法。设x向量的长度取为n，即将积分区间分为n-1段，各段长度为。算出各点的，则矩形法数值积分公式为：2115yx(1,2,,1)ixin(1,2,,1)iyin11niiisyx矩形和梯形定积分公式梯形法的公式为：比较两个公式，它们之间的差别只是。在MATLAB中，把向量中各元素叠加的命令是sum。把向量中各元素按梯形法叠加的命令是trapz。梯形法的几何意义是把被积分的函数的各计算点以直线相联，形成许多窄长梯形条，然后叠加，我们把两种算法都编入同一个程序进行比较。11111222nniiniiiiiyyyyqxyx10.5nyy求面积的数值积分程序exn531fordx=[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),holdoffend程序exn531运行结果程序运行的结果如下：步长dx矩形法解s梯形法解q29108101865815.5841.25816.25.1821.65816.65用解析法求出的精确解为2450/3=816.6666...。dx=2时矩形法和梯形法的积分面积见图5-25.。在曲线的切线斜率为负的情况下，矩形法的积分结果一定偏大，梯形法是由各采样点联线包围的面积，在曲线曲率为负（上凸）时，其积分结果一定偏小，因此精确解在这两者之间。由这结果也能看出，在步长相同时，梯形法的精度比矩形法高。矩形法数字积分的演示程序rsumsMATLAB中有一个矩形法数字积分的演示程序rsums，可以作一个对比。键入rsums('115-x.^2',0,10)就得到右图。图中表示了被积函数的曲线和被步长分割的小区间，并按各区间中点的函数值构成了各个窄矩形面积。用鼠标拖动图下方的滑尺可以改变步长的值，图的上方显示的是这些矩形面积叠加的结果。0246810020406080100115-x.^2:816.99218816MATLAB内的数值定积分函数在实际工作中，用MATLAB中的定积分求面积的函数quad和quadl可以得到比自编程序更高的精度，因为quad函数用的是辛普生法，即把被积函数用二次曲线逼近的算法，而quadl函数采用了更高阶的逼近方法。它们的调用格式如下：Q=QUADL(FUN,A,B,TOL)其中，FUN是表示被积函数的字符串，A是积分下限，B是积分上限。TOL是规定计算的容差，其默认值为1e-6例如，键入S=quad('-x.*x+115',0,10)