高等数学Matlab解法

第5章高等数学问题的MATLAB解法第5章高等数学问题的MATLAB解法5.1函数极限和导数本节习题5.2解析几何和多变量分析本节习题5.3数值积分和微分方程数值解本节习题5.4数列和级数本节习题5.5线性代数本节习题5.6概率论与数理统计本节习题第5章高等数学问题的MATLAB解法5.1函数极限和导数5.1.1单变量函数值的计算和绘图【例5-1-1】要求以0.01s为间隔，求出y的151个点，并求出其导数的值和曲线。解:◆建模可以采取下列两种方法来做:(1)只用主程序编程的方法;(2)编成函数文件，由主程序调用的方法。求导数可采用diff函数对数组y做运算的方法。3πt34sine23y4t第5章高等数学问题的MATLAB解法◆MATLAB程序(1)第一种方法的程序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%求导数并绘制曲线，注意用数值方法求导后，导数数组长度比原函数减少一ylabel(′Dy(t)′)%加标注第5章高等数学问题的MATLAB解法(2)第二种方法的主程序exn511b如下:dt=0.01;t=［0:dt:1.5］;w=4*sqrt(3);y=exn511bf(t，w);Dy=diff(y)/dt;%绘图和加标注的程序略去另要建立一个函数文件exn511bf.m，其内容为:functionxvalues=exn511bf(tvalues，w)%注意编写的函数文件中，其语法对数组w应该能用元素群运算。xvalues=w*exp(-4*tvalues).*sin(w*tvalues+pi/3);第5章高等数学问题的MATLAB解法◆程序运行结果运行这两种程序都得到图5-1所示的曲线。为了节省篇幅，我们没有显示y的数据。以后的各例中还将省略绘图时的标注语句。从本例看，第二种方法似乎更麻烦一些，但它具备模块化的特点。当程序中要反复多次调用此函数，而且输入不同的自变量时，利用函数文件可大大简化编程。我们应该掌握这种方法。两次应用diff函数或用diff(y，2)可以求y的二次导数，读者可自行实践。第5章高等数学问题的MATLAB解法图5-1例5-1-1的曲线第5章高等数学问题的MATLAB解法5.1.2参变方程表示的函数的计算和绘图【例5-1-2】摆线的绘制。如图5-2所示，当圆轮在平面上滚动时，轮上任一点所画出的轨迹称为摆线。如果这一点不在圆周上而在圆内，则生成内摆线;如果该点在圆外，即离圆心距离大于半径，则生成外摆线。对于后一种情况，可由火车轮来想象。其接触轨道的部分，并不是直径最大处，其内侧的直径还要大一些，以防止车轮左右出轨。在这部分边缘上的点就形成外摆线。第5章高等数学问题的MATLAB解法图5-2摆线的生成第5章高等数学问题的MATLAB解法解：◆建模概括几种情况，其普遍方程可表示为:xA=rt-RsintyA=r-Rcost可由这组以t为参数的方程分析其轨迹。第5章高等数学问题的MATLAB解法◆MATLAB程序t=0:0.1:10;%设定参数数组r=input(′r=′)，R=input(′R=′)%输入常数x=r*t-R*sin(t);%计算x，yy=r-R*cos(t);plot(x，y)，axis(′equal′)%绘图第5章高等数学问题的MATLAB解法◆程序运行结果设r=1，令R=r，R=0.7及R=1.5时得到的摆线、内摆线和外摆线都绘于图5-3中。为了显示摆线的正确形状，x、y坐标保持等比例是很重要的，因此程序中要加axis(′equal′)语句。第5章高等数学问题的MATLAB解法图5-3摆线的绘制第5章高等数学问题的MATLAB解法5.1.3曲线族的绘制【例5-1-3】三次曲线的方程为y=ax3+cx，试探讨参数a和c对其图形的影响。解:◆建模因为函数比较简单，因此可以直接写入绘图语句中，用循环语句来改变参数。注意坐标的设定方法，以得到适于观察的图形。给出的程序不是唯一的，例如也可用fplot函数等。读者可自行探索其他编法。第5章高等数学问题的MATLAB解法◆MATLAB程序x=-2:0.1:2;%给定x数组，确定范围及取点密度subplot(1，2，1)%分两个画面绘图forc=-3:3%取不同的c循环plot(x，x.^3+c*x)，holdon，end，gridsubplot(1，2，2)fora=-3:3%取不同的a循环plot(x，a*x.^3+x)，holdon，end，grid，holdoff用直接在图形窗内编辑的方法可在图内标注字符。第5章高等数学问题的MATLAB解法◆程序运行结果程序运行结果见图5-4，其中a和c均从-3取到3，步长为1。第5章高等数学问题的MATLAB解法图5-4c和a取不同值时y=ax3+cx(a)a=1，c取不同值;(b)c=1，a取不同值第5章高等数学问题的MATLAB解法5.1.4极限判别用MATLAB来表达推理过程是比较困难的，因为它必须与实际的数值联系起来。比如无法用无穷小和高阶无穷小的概念，只能用e-10、e-20等数值。不过极限的定义恰恰用了δ和ε这些数的概念，因此不难用程序表达。【例5-1-4】极限的定义和判别。解:◆建模对于函数y=f(x)，当任意给定一个正数ε时，有一个对应的正数δ0＜|xc-x|δ时,|A-f(x)|＜ε其中，A是f(x)在x→xc时的极限，如果找不到这样的δ，A就不是它的极限。只考虑左极限时，因为xc-x必为正数，所以可去掉绝对值符号。第5章高等数学问题的MATLAB解法◆MATLAB程序检验极限是否正确的程序disp(′A是否是f(xc)的极限?′)fxc=input(′f(x)的表达式为，例如sin(x)/x′，′s′)，%输入函数表达式A=input(′A=，例如A=1′)，%输入极限值xc=input(′xc=，例如xc=0′)，%输入对应的自变量值flag=1;delta=1;x=xc-delta;n=1;%初始化whileflag==1epsilon=input(′任给一个小的数ε=′)%任意给出εwhileabs(A-eval(fxc))epsilondelta=delta/2;x=xc-delta;%找δifabs(delta)epsdisp(′找不到δ′)，n=0;break%找不到δ，跳出内循环end，end第5章高等数学问题的MATLAB解法ifn==0disp(′极限不正确′)，break，end，%极限不正确，跳出外循环delta，disp(′极限正确′)%找到了δflag=input(′再试一个ε吗?再试按1，不试按0或任意数字键′)，%再试end第5章高等数学问题的MATLAB解法◆程序运行结果试判断:(1)f(x)=x2-8在x→xc=3时是否以1.001为左极限？(2)f(x)=sin(x)/x在x→xc=0时是否以1为左极限？对于(1)，将得出“极限不正确”的结果；而对于(2)，将得出“极限可能正确”的结果。读者可思考为什么要加“可能”二字而不能给出肯定的答复。第5章高等数学问题的MATLAB解法5.1.5画出曲线并求左右极限【例5-1-5】画出的图形，找到它的间断点x0，判断它在间断点附近的取值，并求它在x=0处的值或极限值。xsinxcosx1f(x)第5章高等数学问题的MATLAB解法解:◆程序及运行结果用数值计算方法来求。先绘图，程序如下：ezplot(′(1-cos(x))./x./sin(x)′)gridon运行程序得出图5-5。从图上可见，此函数在x0=±π处有间断点。在间断点左右，函数f(x)分别趋向于±∞。现在尝试用数组方法计算以下三个特征点的值，令x1=［-pi，0，pi］y1=(1-cos(x1))./x1./sin(x1)得到警告:Warning:Dividebyzero及y1=1.0e+015*［5.1984，NaN，5.1984］第5章高等数学问题的MATLAB解法图5-5例5-1-5的曲线图形第5章高等数学问题的MATLAB解法可见左右x0=±π处的两个数都是很大的，而中间x=0处则是非数(NaN)。这是由于在该点出现了函数y=0/0所造成的。从理论上说，应该通过罗必达法则，分别对分子和分母求导，再求其极限。但求导也有点麻烦，在用计算机解题时，宁可用同一个程序多算几个点。为此，对这三个特殊点的左右邻域各取一点，都做计算，共计算六个点。令e1=0.001，x2=[-pi-e1，-pi+e1，0-e1，0+e1，pi-e1，pi+e1]y2=(1-cos(x2))./x2./sin(x2)第5章高等数学问题的MATLAB解法得到的结果是y2=-636.4171636.82240.50000.5000636.8224-636.4171这就证实了我们的判断大体上是对的，即在x0=±π的间断点两侧，函数f(x)分别趋向于-∞和∞，而在x=0处，f(x)趋向于0.5。第5章高等数学问题的MATLAB解法◆讨论用数值方法虽然方便，但邻域大小的选择有时会影响计算结果。选得大了，结果的偏差比较大，比如这里算出的636离无穷大还差得远，但也不能选得太小。几个很小的数之间作除法所造成的误差往往极其可观，比如在这道题中把邻域选为最小数eps，就得不到正确的结果。一般来说，有了全局的曲线，从图上就可判断大致的结果，因为ezplot函数的自变量增量是上下限之间分成300份，如果使间隔值缩小，且在临界点附近函数值看不出大的起伏，那这个结果就可以相信了。第5章高等数学问题的MATLAB解法如果求x在无穷大处的极限，可以进行自变量替换，用x1=1/x把f(x)变为f1(x1),然后研究以x1为自变量的函数的曲线，并研究x1=0附近的极限和特性。比如，若f(x)=xex，可用下面的方法求。令x1=1/x,得到，写出以下程序:x1=linspace(-pi，pi，1001);plot(x1，exp(1./x1)./x1)，gridonaxis(［-pi，pi，-2，10］)f(x)limx)(xfxe)f(1/xy1111/x11第5章高等数学问题的MATLAB解法运行后得出图5-6，然后求x1=0的邻域y值，键入:x1=［-eps，0，eps］;y=exp(1./x1)./x1Warning:Dividebyzero.y=0InfInf这说明x1取-eps即x趋向-Inf时，y的极限为零;而在x趋向Inf时，y趋向Inf。用MATLAB中符号数学的语句可以直接求解如下：symsx，y=limit(x*exp(x)，x，Inf)得到y=Inf,symsx，y=limit(x*exp(x)，x，-Inf)得到y=0。第5章高等数学问题的MATLAB解法图5-6函数y=f1(x1)的图形第5章高等数学问题的MATLAB解法【例5-1-6】(a)画出g(x)=xsin(1/x)，估计;(b)画出k(x)=sin(1/x)，估计;(c)比较g和k在原点附近的特性，有什么相同点和不同点？解:先画出这两个函数的概略曲线形状，程序为subplot(1，2，1)，ezplot(′x.*sin(1./x)′)，subplot(1，2，2)，ezplot(′sin(1./x)′)得到图5-7，其中g(x)对应于图5-7(a)，从中大体可以看出，它在原点附近摆动的幅度似乎在减小;k(x)对应于图5-7(b)，该曲线是反对称的，在原点附近似乎有大幅度的摆动，需要把图形进一步放大。g(x)limxk(x)limx第5章高等数学问题