第3章微积分学

1、第3章 Matlab在微积分学中的应用在这一章里，我们将利用MATLAB解决高等数学中的许多计算问题。我们还将学会建立一般的函数在一定程度上解决理论上的计算问题。3.1 极限数学中的极限问题类型大概有：极限数学形式的例子：数列的极限有：一元函数的极限有：二元函数的极限有：以上是数学中出现的不同形式的极限。但求极限的命令在MATLAB中是相同的，都由命令limit来完成。其常用形式有：limit(F,x,a) 计算limit(F,a) 符号表达式F中由命令findsym(F)返回独立的变量v，计算limit(F) 符号确定同上，设为v, 计算limit(F,x,a,right) 计算limit(

2、F,x,a,left) 计算注意：求一极限的基本步骤为：1.把数学形式的表达式转化为Matlab的表达式。当然表达式由多个符号组成，并且符号都不是代表某一具体的数值，其值是一个变动的量，即所谓的符号变量。为此，我们要建立一些必要的符号对象，如：符号变量，符号表达式等。2.使用求极限的命令limit来对符号表达式进行计算，其中要指明是对哪一个符号变量做极限运算；若不指定，则符号变量由命令findsym(F)来指定。3.要指定符号变量的“极限点”以及趋于“极限点”的方向；若不指定“极限点”的值，则缺省地认为是0点。4.总之一点，极限作的都是符号运算！例3-1求解：>> syms n m

3、 xf=(6*n2-n+1)/(n3-n2+2)g=(1+m*x)n-(1+n*x)m)/x2h=(sqrt(1+x)-2)/(x-3)>> lim_f=limit(f,n,inf)>> lim_g=limit(g,x,0) %或lim_g=limit(g)>> lim_h=limit(h,x,3,'right')运算结果是： lim_f =0 lim_g =-1/2*m2*n+1/2*n2*m lim_h =1/4例3-2 求由于求二重极限在数学上没有方法可循，因此在Matlab中还没有一个统一的命令能求一个一般的二重极限，只能求在理论上已

4、经证明与路径无关的函数，即把二重极限化成二次极限来计算：解：>> syms x y>> f_xy='log(x+exp(y)/sqrt(x2+y2)'>> lim_f_xy=limit(limit(f_xy,x,1),y,0)运行的结果是：lim_f_xy =log(2)3.2微积分高等数学中的微积分是大学数学中的重要部分，内容庞杂，应用广泛，归纳起来有：函数的导数，函数的积分，泰勒展式，-函数，欧拉函数等。下面先讲函数的导数。3.2.1函数的导数数学中的导数类型大概有：导数例如：一元函数的导数：y=esinx-7cosx+5x,y=二元函数

5、的导数：z=arctan,求，z=ln(复合函数的导数：y=f(x)=a+bx+f(x)= 隐函数的导数：arctgF(x,x+y,x+y+z)=0,求。数学上虽然有导数与偏导数之分，但它们在Matlab中有统一的命令diff。其常使用形式为：diff(f,v)计算diff(f,v,n)计算diff(f)用命令findsym找出表达式f中的独立变量，设为v,即v=findsym(F),计算diff(f,n)独立变量确定同上，设为v,计算例3-3 已知y=计算解：>> x=sym('x')>> y=(sqrt(x)+cos(x)/(x-1)-7*x2;&g

6、t;> Dy=diff(y) >>D3y=diff(y,3)运行结果为：Dy =(1/2/x(1/2)-sin(x)/(x-1)-(x(1/2)+cos(x)/(x-1)2-14*xD3y =(3/8/x(5/2)+sin(x)/(x-1)-3*(-1/4/x(3/2)-cos(x)/(x-1)2+6*(1/2/x(1/2)-sin(x)/(x-1)3-6*(x(1/2)+cos(x)/(x-1)4例3-4 已知z=ln(。解：>> x=sym('x');y=sym('y');z=log(sqrt(x)+sqrt(y)；result

7、=x*diff(z,x)+y*diff(z,y)；simple (result)计算的结果为：result=1/2例3-5已知y=ln(x),x=t。解：在复合函数的计算中，一定要注意变量赋值的先后顺序>> syms t x>> x=t2*sin(t);>> y=log(x3);>> dydt=diff(y,t)>> d2yd2t=diff(dydt,t)/diff(x,t)计算结果为： dydt =(6*t5*sin(t)3+3*t6*sin(t)2*cos(t)/t6/sin(t)3d2yd2t =(30*t4*sin(t)3+3

8、6*t5*sin(t)2*cos(t)+6*t6*sin(t)*cos(t)2-3*t6*sin(t)3)/t6/sin(t)3-6*(6*t5*sin(t)3+3*t6*sin(t)2*cos(t)/t7/sin(t)3-3*(6*t5*sin(t)3+3*t6*sin(t)2*cos(t)/t6/sin(t)4*cos(t)/(2*t*sin(t)+t2*cos(t)例3-6 已知arctg解：设方程F（x,y,z）=0确定了函数z=z(x,y),则=-，再用类似的计算方法，可以计算。>> syms x y z>> F=atan(y+z)/x)-log(x+y+z)

9、>> dFdx=diff(F,x)>> dFdz=diff(F,z)>> dZdx=-dFdx/dFdz>> dZdxdy=diff(dZdx,y)>> dZdxdz=diff(dZdx,z)>> d2Zdxdy=-dZdxdy/dZdxdz计算结果为：dFdx =-(y+z)/x2/(1+(y+z)2/x2)-1/(x+y+z)dFdz =1/x/(1+(y+z)2/x2)-1/(x+y+z)dZdx =(y+z)/x2/(1+(y+z)2/x2)+1/(x+y+z)/(1/x/(1+(y+z)2/x2)-1/(x+y+

10、z)dZdxdy =(1/x2/(1+(y+z)2/x2)-2*(y+z)2/x4/(1+(y+z)2/x2)2-1/(x+y+z)2)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)-(y+z)/x2/(1+(y+z)2/x2)+1/(x+y+z)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)2*(-2/x3/(1+(y+z)2/x2)2*(y+z)+1/(x+y+z)2)dZdxdz =(1/x2/(1+(y+z)2/x2)-2*(y+z)2/x4/(1+(y+z)2/x2)2-1/(x+y+z)2)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)-(y+z

11、)/x2/(1+(y+z)2/x2)+1/(x+y+z)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)2*(-2/x3/(1+(y+z)2/x2)2*(y+z)+1/(x+y+z)2)d2Zdxdy =(-(1/x2/(1+(y+z)2/x2)-2*(y+z)2/x4/(1+(y+z)2/x2)2-1/(x+y+z)2)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)+(y+z)/x2/(1+(y+z)2/x2)+1/(x+y+z)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)2*(-2/x3/(1+(y+z)2/x2)2*(y+z)+1/(x+y+z)2)/

12、(1/x2/(1+(y+z)2/x2)-2*(y+z)2/x4/(1+(y+z)2/x2)2-1/(x+y+z)2)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)-(y+z)/x2/(1+(y+z)2/x2)+1/(x+y+z)/(1/x/(1+(y+z)2/x2)-1/(x+y+z)2*(-2/x3/(1+(y+z)2/x2)2*(y+z)+1/(x+y+z)2)3.2.2 一重积分3.2.2.1 不定积分不定积分是高等数学中的基本运算。数学中的许多不定积分的原函数不能用初等函数来表示，因此在Matlab中，有些函数的原函数是求不出来的。Matlab提供了几个用于求函数不定积分的

13、命令，如：int，intwave，fnint (用于专门函数的积分)。其中函数int是最常用的，功能强大的积分函数。下面我们介绍函数int的使用格式：F = int (f)F = int (f, var)参数说明：·f：被积函数；·var：积分变元；若没有指定，则对变量v = findsym (f)积分；·F：函数f的原函数。F中没有常数C。例3-7 计算 ，解：syms x y z>> f=exp(x*y+z);>> F1=int(f) F1 = 1/y*exp(x*y+z)>> F2=int(f,z) F2 = exp(x*

14、y+z)3.2.2.2 定积分及广义积分数学中的定积分有时和不定积分使用的方法不同。我们这里所讲的广义积分其积分区间端点中至少有一个为无穷，而函数在该区间内是连续的。在Matlab中，只需在命令int中加入积分上下限即可。即F = int(f,a,b)F = int(f,var,a,b)参数说明：·f：被积函数；·a,b：积分上下限；a, b可为无穷大（若为负无穷大，则为：）；·F：函数的积分值。有时为无穷大。定积分的计算，也可调用函数quad ( )或quad8 ( )，其格式为：y, n = quad (f , a, b, tol)y, n = quad8 (

15、f , a, b, tol)参数说明：·f：自定义的被积函数；·a,b：积分上下限；·n：计算过程中调用被积函数f的次数；tol：精度要求，quad ( )的默认值为：tol = 1e-3，quad8 ( )的默认值为：tol = 1e-6；·quad ( )函数采用的是Simpson方法计算定积分的近似值；·quad8 ( )函数采用的是Newton Cotes 方法计算定积分的近似值，其精度比前者更高。例3-8 计算， 解：f=1/(x2+2*x+3);>> F1=int(f,2,pi) F1 = 1/2*atan(1/2*2(

16、1/2)*(pi+1)*2(1/2)-1/2*atan(3/2*2(1/2)*2(1/2) >> F2=int(f,-inf,inf) F2 = 1/2*pi*2(1/2)例3-9 计算解：先定义函数，文件名：f.mfunction y=f(x)y=1/sqrt(2*pi)*exp(-x.2/2);保存后，在命令窗口键入format long>> y,n=quad('f',-3,3)则显示结果为：y = 0.99729991863154n = 57 （表示被积函数f的调用次数）3.2.3 多重积分3.2.3.1 二重积分一种方法：由于积分区域的多样性和复

17、杂性，通过人为地结合积分区域的图形，把二重积分转化为二次积分，用int命令进行计算；另一种方法：利用Matlab提供的函数dblquad和quad2dggen计算，其调用格式如下：·二重数值积分（在矩形区域上），使用dblquad函数，格式为： y = dblquad (f, x_m, x_M, y_m, y_M, tol)·一般二重积分（在数值积分工具箱（NIT，即Numerical Integration Toolbox）中提供了quad2dggen函数，该函数可从网上下载，网址：contrib./v4/inte-gration/nit.zip，安装之后，将其路径加载到

18、Matlab路径下），格式为： I = quad2dggen（被积函数名，下限函数名，上限函数名，y_m，y_M，tol）例3-10 计算其中D为直线y = x和抛物线y = x2所围部分。解：由数学方法可得：（1）画出积分区域示意图>>syms x y>> f=(2-x-y)/2;>> y1=x;>> y2=x2;>> ezplot(y1);hold on>> ezplot(y2);>> axis(0,2,0,2)（2）确定积分限a=fzero('x-x2',0)a = 0>> b

19、=fzero('x-x2',1)b = 1这是第2次积分的上下限。（3）积分运算f_dy=int(f,y,x2,x) %先对y积分f_dy = x-5/4*x2-1/2*x*(x-x2)+1/4*x4>> I=int(f_dy,a,b) %再对x积分I = 11/120例3-11 试求下面的二次积分解：先定义函数，函数名为：my2dfun.mfunction z=my2dfun(x,y)global kk;kk=kk+1; %定义全局变量，测定被积函数调用次数z=exp(-x.2/2).*sin(x.2+y);保存后在命令窗口执行下列命令：clear>>

20、 global kk;kk=0;>> y=dblquad('my2dfun',-2,2,-1,1),kky =1.57449318974494kk = 1038 （表示被积函数调用次数）3.2.3.2 三重积分三重积分的过程与二重积分过程相同，也是把三重积分转化成三次积分来计算，只不过在确定积分限时更繁琐而已。在确定积分限时一般要结合三维的积分区域，所以一般先画出积分区域。例3-12 计算其中区域V为三个坐标面及平面：所围成的闭区域。解：由数学方法可得： （1）画出积分区域示意图。由于图形比较简单，所以我们没有画出。（并且Matlab没有画三维符号函数的命令）sym

21、s x y z f=x;（2）确定积分限。a=0;b=1;phi1=0;phi2=(1-x)/2;psi1=0;psi2=1-x-2*y;（3）积分计算。f_dz=int(f,z,psi1,psi2)f_dzdy=int(f_dz,y,phi1,phi2)I=int(f_dzdy,x,a,b)计算结果为：f_dz = x*(1-x-2*y)f_dzdy =-x*(1/2-1/2*x)2+x*(1-x)*(1/2-1/2*x)I =1/483.2.4 Taylor展式数学定理表明：若f (x)在x = 0点附近有直到n+1阶连续导数，那么：其中： 实际上是要将函数f (x)表示成xn（n从0到无

22、穷）的和的形式。这完全可以用Matlab提供的命令taylor来完成展开工作。其常用的使用形式为：·taylor (f)·taylor (f, n)·taylor (f, v)·taylor (f,n, a)参数说明：f：待展开的函数表达式，可以不用单引号生成；n：把函数展开到n阶，若不包含n，则缺省地展开到6阶；v：对函数f中的变量v展开，若不包含v，则对变量v = findsym (f)展开；a：Taylor展式的扩充功能，对函数f在x = a点展开。例3-13 计算（1）把y = e x展开到6阶；（2）把y = lnx在x = 1点展开到6阶；（

23、3）把y = sinx在x =/2点展开到6阶；（4）把y = x t关于变量t展开到3阶。解：syms x ty1=taylor(exp(-x)y2=taylor(log(x),6,1)y3=taylor(sin(x),pi/2,6)y4=taylor(xt,3,t)运行结果为：y1 = 1-x+1/2*x2-1/6*x3+1/24*x4-1/120*x5y2 = x-1-1/2*(x-1)2+1/3*(x-1)3-1/4*(x-1)4+1/5*(x-1)5y3 = 1-1/2*(x-1/2*pi)2+1/24*(x-1/2*pi)4y4 = 1+log(x)*t+1/2*log(x)2*t

24、23.2.5 Fourier级数数学定理表明：设函数f (x)已经展开为全区间上的一致收敛的三角级数：其中： 虽然Matlab中还没有提供一个专门的命令用于求解Fourier级数，但根据上面的数学公式，也可以编写一个函数文件sfour.m，用于求函数f在区间上的Fourier级数的系数：function a0,ak,bk=sfour(f)syms x na0=int(f,-pi,pi)/pian=int(f*cos(n*x),-pi,pi)/pi;an=simple(an)bn=int(f*sin(n*x),-pi,pi)/pi;bn=simple(bn)有了上面的函数，我们可以求得一些函数的

25、Fourier级数的系数。例3-14 求函数y = x的Fourier级数的系数解：调用上面定义的函数sfour来计算：>>syms x>> f=x;>> a0,ak,bk=sfour(f)a0 =0an =0bn =2/n2/pi*sin(pi*n)-2/n*cos(pi*n)注：在Matlab中，化简命令不能把sin(n*pi)和cos(n*pi)转化为(-1)n。所以bn的结果较长。3.2.6 梯度设u = u(x, y, z)是一数量函数，点p (x p, y p, z p)是等高面u = u (x, y, z) = C上的任一点，则称向量为函数u

26、= u (x, y, z)在p点的梯度。1. 近似梯度函数gradient。使用形式： （1）FX, FY = gradient (F) （2）FX, FY = gradient (F, H) （3）FX, FY = gradient (F, HX, HY) （4）FX, FY, FZ = gradient (F) （5）FX, FY, FZ = gradient (F, HX, HY, HZ) （6）FX, FY, FZ, = gradient (F, )参数说明： （1）F：待求梯度的数值矩阵。F可以为向量、二维矩阵、三维矩阵； （2）FX，FY，FZ：矩阵F在x、y、z方向上的数值梯度；

27、 （3）H：H为一标量，作为各个方向上各点之间的步长； （4）HX，HY，HZ：矩阵F在x、y、z方向上的具体步长。HX，HY，HZ可为标量或与矩阵F各个方向上同维的向量。例3-15 计算函数数值梯度，且以图形显示。解：x,y=meshgrid(-2:.2:2,-2:.2:2); %由两个向量生成的矩阵>> z=x.*exp(-x.2-y.2);>> px,py=gradient(z,0.2,0.2);>> contour(z) %等高线图>> hold on>> quiver (px,py) %向量图运行结果为：数值梯度图2. 函数

28、梯度和方向导数jacobian。使用形式：jacobian (f, v) 参数说明： f：函数向量或标量，当f为标量时，jacobian (f, v) = gradient (f); v：自变量向量或者单个变量。例3-16求：（1）u = x y z在点播M (1, 1, 1)处的梯度；（2）u = x2+2y2+3z2+xy+3x-2y-6z在点O (0, 0, 0)及A (1, 2, 3)处的梯度大小。syms x y z >> u1=x*y*z;>> u2=x2+2*y2+3*z2+x*y+3*x-2*y-6*z;>> v=x,y,z;>>

29、 J1=jacobian(u1,v);>> J2=jacobian(u2,v);J1_M=subs(subs(subs(J1,x,1),y,1),z,1)%subs(x, 1) 用1替换x，计算其值J1_M = 1 1 1>> J2_O=subs(subs(subs(J2,x,0),y,0),z,0)J2_O = 3 -2 -6>> J2_A=subs(subs(subs(J2,x,1),y,2),z,3)J2_A = 7 7 123.3 最小二乘法在生产实践中，常常需要根据实际测量得到的一系列数据找出函数关系，通常叫做配曲线或找经验公式。本节我们介绍一种找

30、直线的方法，它是广泛使用的一种处理数据的方法。例3-17 在某实验中测得数据如下：X104180190177147134150191204121Y100200210185155135170205235125由此要推出x和y的函数关系：y = f (x)。解：先把这些数据点描绘出来，如图所示，观察x和y大概满足的函数关系。先把x和y进行适当的调整，使自变量x的值从小到大排列。>>x=104 180 190 177 147 134 150 191 204 121;>> y=100 200 210 185 155 135 170 205 235 125;>> x,

31、i=sort(x)x = 104 121 134 147 150 177 180 190 191 204i = 1 10 6 5 7 4 2 3 8 9>> y=y(i)y = 100 125 135 155 170 185 200 210 205 235>> plot(x,y,'r*')>> hold on这些数据点大致分布在一条直线上，即x和y有线性函数关系：y = ax+b，由数学推断过程得：， 其中n为x和y的长度。用Matlab计算为：x=104 180 190 177 147 134 150 191 204 121;>>

32、; y=100 200 210 185 155 135 170 205 235 125;>> x,i=sort(x) %以升序排序，i为元素在原向量中的位置。x = 104 121 134 147 150 177 180 190 191 204i = 1 10 6 5 7 4 2 3 8 9>> y=y(i)y = 100 125 135 155 170 185 200 210 205 235>> plot(x,y,'r*')>> hold on>> n=length(x);>> a_den=n*sum(x

33、.*y)-sum(x)*sum(y);>> b_den=sum(y).*sum(x.2)-sum(x).*sum(x.*y);>> ab_num=n*sum(x.2)-(sum(x)2;>> a=a_den/ab_numa = 1.2673>> b=b_den/ab_numb = -30.5143>> x=100:0.5:220;>> y=a*x+bplot(x,y)运算结果为：最小二乘法结果3.4 微积分的应用 例3-18 求通过点(1, 1)的曲线方程，使其在区间1，x上对应的曲边梯形面积等于该曲线上点M (x, y)

34、的两个坐标之积的两倍（x>0, y>0）。解：设曲线为y = y (x)，按所给条件，可得对上式两边求导，得微分方程y = 2 (y-xy)，或y= y(2x)，且y (1) = 1。对此，可用dsolve函数求解。syms x yy=dsolve('Dy+y/(2*x)=0','y(1)=1','x')计算结果为：y =1/x(1/2)或syms t y>> y=dsolve('Dy+y/(2*t)=0','y(1)=1')计算结果为：y =1/t(1/2)3.5 方程（组）的求解数学计算

35、中有很大一部分是求解方程（组），有些方程（组）可以得到精确的解析解，而有些却很难做到，只能求得满足一定精度的数值解。对于线性方程组的解，我们在后面章节中讨论。这一节主要是讨论一般的代数方程（组）的解。在Matlab中解方程组的命令有linsolve和solve，其中命令linsolve专门用于求解线性方程组，而命令solve可适用于所有的代数方程（组）。其常用形式有：1. x = linsolve (A, B) ss = solve (s) ss = solve (s, v)参数说明：·s：包含方程（一个）等式的字符串（可以是函数名，或者是描述方程的字符串）；·v：方程s中

36、的一个变量；·ss：若是第一种情形，对s中默认的变量求解，结果赋给ss；若是第二种情形，对s中指定的变量v求解，结果赋给ss；·x：x是矩阵方程Ax = B的解x = sym(A)sym(B) = A-1*B。例3-19 求解：（1）Ax = B，其中A，B； （2）psin(2x+t) = q，其中t为未知参数。解：A=2 5;1 3;B=4 -6;2 1;>> x=linsolve(A,B)x = 2, -23 0, 8 >> solve('p*sin(2*x+t)=q','t') ans =-2*x+asin(q

37、/p)2. ss = solve (s1, s2, , sn) ss = solve (s1, s2, , sn, v1, v2, , vm) x1, x2, , xn = solve (s1, s2, , sn) x1, x2, , xn = solve (s1, s2, , sn, v1, v2, , vm)参数说明：·s1, s2, , sn：包含方程组（n个）等式的字符串（可以是函数名或者是描述方程组的字符串表达式）·x1, x2, , xn：若是第一种情形，对n个方程s1, s2, , sn中默认的变量求解；若是第二种情形，对n个方程s1, s2, , sn中指定

38、的变量求解；结果分别赋给向量x1, x2, , xn中的分量；·v1, v2, , vm：方程组s1, s2, , sn中的m个变量；·ss：在第一种情形中，对n个方程s1, s2, , sn中默认的变量求解；若是第二种情形，对n个方程s1, s2, , sn中指定的变量求解；n与m不一定相等，即可能是超定方程组或是不定方程组，结果赋给出ss，ss为一个解向量。例3-20 求下列方程组的解（1） （2）（3） （4）解：x1,y1=solve('x2+x*y+y=3','x2-4*x+3=0')x2,y2=solve('a*x2+y2

39、=0','x-y=1')x3,y3,z3=solve('x*y2+z2=0','y-z=1','x2-5*x+6=0')x,y=solve('sin(x+y)-exp(x)*y=0','x2-y=2')计算结果为：x1 = 1 3y1 = 1 -3/2x2 = 1/2/(a+1)*(-2*a+2*(-a)(1/2)+1 1/2/(a+1)*(-2*a-2*(-a)(1/2)+1y2 = 1/2/(a+1)*(-2*a+2*(-a)(1/2) 1/2/(a+1)*(-2*a-2*(-a)(1/

40、2)x3 = 2 2 3 3y3 = 1/3+1/3*i*2(1/2) 1/3-1/3*i*2(1/2) 1/4+1/4*i*3(1/2) 1/4-1/4*i*3(1/2)z3 = -2/3+1/3*i*2(1/2) -2/3-1/3*i*2(1/2) -3/4+1/4*i*3(1/2) -3/4-1/4*i*3(1/2)x =y = 34.208227234306296508646214438330由于第四个方程组没有精确的解析解，所以给出的是数值形式的解。3.6 函数极值在高等数学中，求函数的极值（最值）是函数导数的应用之一。不过高等数学中的求值过程比较麻烦，而在Matlab中却显得轻松自如。因为Matlab有两个功能强大的求值函数fminbnd、fminsearch。3.6.1 一元函数的极值函数fminbnd是用于专门求单变量函数的最小值。使用形式为：x = fminbnd(fun,x1,x2)x = fminbnd(fun,x1,x2,options)x = fminbnd(fun,x1,x2,options,p1,p2,.)x,fval = fminbnd(.)x,fval,exitflag = fminbnd(.)x,fval,exitflag,output = fminbnd(.)参数说明：·fun：目标