第3章微积分学.docx
《第3章微积分学.docx》由会员分享，可在线阅读，更多相关《第3章微积分学.docx（27页珍藏版）》请在冰豆网上搜索。
第3章微积分学
第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（F,x,a,’left’）计算
注意：
求一极限的基本步骤为：
1.把数学形式的表达式转化为Matlab的表达式。
当然表达式由多个符号组成，并且符号都不是代表某一具体的数值，其值是一个变动的量，即所谓的符号变量。
为此，我们要建立一些必要的符号对象，如：
符号变量，符号表达式等。
2.使用求极限的命令limit来对符号表达式进行计算，其中要指明是对哪一个符号变量做极限运算；若不指定，则符号变量由命令findsym（F）来指定。
3.要指定符号变量的“极限点”以及趋于“极限点”的方向；若不指定“极限点”的值，则缺省地认为是0点。
4.总之一点，极限作的都是符号运算！
例3-1求
解：
>>symsnmx
f=（6*n^2-n+1）/（n^3-n^2+2）
g=（（1+m*x）^n-（1+n*x）^m）/x^2
h=（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*m^2*n+1/2*n^2*m
lim_h=
1/4
例3-2求
由于求二重极限在数学上没有方法可循，因此在Matlab中还没有一个统一的命令能求一个一般的二重极限，只能求在理论上已经证明与路径无关的函数，即把二重极限化成二次极限来计算：
解：
>>symsxy
>>f_xy='log（x+exp（y））/sqrt（x^2+y^2）';
>>lim_f_xy=limit（limit（f_xy,x,1）,y,0）
运行的结果是：
lim_f_xy=
log
（2）
3.2微积分
高等数学中的微积分是大学数学中的重要部分，内容庞杂，应用广泛，归纳起来有：
函数的导数，函数的积分，泰勒展式，
-函数，欧拉函数等。
下面先讲函数的导数。
3.2.1函数的导数
数学中的导数类型大概有：
导数
例如：
一元函数的导数：
y=e
sinx-7cosx+5x
y=
二元函数的导数：
z=arctan
求
，
z=ln（
复合函数的导数：
y=f（x）=a+bx+
f（x）=
隐函数的导数：
arctg
F（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*x^2;
>>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*x
D3y=
（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=x*diff（z,x）+y*diff（z,y）；
simple（result）
计算的结果为：
result=
1/2
例3-5已知y=ln（x
）,x=t
。
解：
在复合函数的计算中，一定要注意变量赋值的先后顺序
>>symstx
>>x=t^2*sin（t）;
>>y=log（x^3）;
>>dydt=diff（y,t）
>>d2yd2t=diff（dydt,t）/diff（x,t）
计算结果为：
dydt=
（6*t^5*sin（t）^3+3*t^6*sin（t）^2*cos（t））/t^6/sin（t）^3
d2yd2t=
（（30*t^4*sin（t）^3+36*t^5*sin（t）^2*cos（t）+6*t^6*sin（t）*cos（t）^2-3*t^6*sin（t）^3）/t^6/sin（t）^3-6*（6*t^5*sin（t）^3+3*t^6*sin（t）^2*cos（t））/t^7/sin（t）^3-3*（6*t^5*sin（t）^3+3*t^6*sin（t）^2*cos（t））/t^6/sin（t）^4*cos（t））/（2*t*sin（t）+t^2*cos（t））
例3-6已知
arctg
解：
设方程F（x,y,z）=0确定了函数z=z（x,y）,则
=-
，再用类似的计算方法，
可以计算
。
>>symsxyz
>>F=atan（（y+z）/x）-log（x+y+z）
>>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）/x^2/（1+（y+z）^2/x^2）-1/（x+y+z）
dFdz=
1/x/（1+（y+z）^2/x^2）-1/（x+y+z）
dZdx=
（（y+z）/x^2/（1+（y+z）^2/x^2）+1/（x+y+z））/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））
dZdxdy=
（1/x^2/（1+（y+z）^2/x^2）-2*（y+z）^2/x^4/（1+（y+z）^2/x^2）^2-1/（x+y+z）^2）/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））-（（y+z）/x^2/（1+（y+z）^2/x^2）+1/（x+y+z））/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））^2*（-2/x^3/（1+（y+z）^2/x^2）^2*（y+z）+1/（x+y+z）^2）
dZdxdz=
（1/x^2/（1+（y+z）^2/x^2）-2*（y+z）^2/x^4/（1+（y+z）^2/x^2）^2-1/（x+y+z）^2）/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））-（（y+z）/x^2/（1+（y+z）^2/x^2）+1/（x+y+z））/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））^2*（-2/x^3/（1+（y+z）^2/x^2）^2*（y+z）+1/（x+y+z）^2）
d2Zdxdy=
（-（1/x^2/（1+（y+z）^2/x^2）-2*（y+z）^2/x^4/（1+（y+z）^2/x^2）^2-1/（x+y+z）^2）/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））+（（y+z）/x^2/（1+（y+z）^2/x^2）+1/（x+y+z））/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））^2*（-2/x^3/（1+（y+z）^2/x^2）^2*（y+z）+1/（x+y+z）^2））/（（1/x^2/（1+（y+z）^2/x^2）-2*（y+z）^2/x^4/（1+（y+z）^2/x^2）^2-1/（x+y+z）^2）/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））-（（y+z）/x^2/（1+（y+z）^2/x^2）+1/（x+y+z））/（1/x/（1+（y+z）^2/x^2）-1/（x+y+z））^2*（-2/x^3/（1+（y+z）^2/x^2）^2*（y+z）+1/（x+y+z）^2））
3.2.2一重积分
3.2.2.1不定积分
不定积分是高等数学中的基本运算。
数学中的许多不定积分的原函数不能用初等函数来表示，因此在Matlab中，有些函数的原函数是求不出来的。
Matlab提供了几个用于求函数不定积分的命令，如：
int，intwave，fnint（用于专门函数的积分）。
其中函数int是最常用的，功能强大的积分函数。
下面我们介绍函数int的使用格式：
F=int（f）
F=int（f,var）
参数说明：
·f：
被积函数；
·var：
积分变元；若没有指定，则对变量v=findsym（f）积分；
·F：
函数f的原函数。
F中没有常数C。
例3-7计算
，
解：
symsxyz
>>f=exp（x*y+z）;
>>F1=int（f）
F1=
1/y*exp（x*y+z）
>>F2=int（f,z）
F2=
exp（x*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（‘f’,a,b,tol）
参数说明：
·f：
自定义的被积函数；
·a,b：
积分上下限；
·n：
计算过程中调用被积函数f的次数；tol：
精度要求，quad（）的默认值为：
tol=1e-3，quad8（）的默认值为：
tol=1e-6；
·quad（）函数采用的是Simpson方法计算定积分的近似值；
·quad8（）函数采用的是NewtonCotes方法计算定积分的近似值，其精度比前者更高。
例3-8计算
，
解：
f=1/（x^2+2*x+3）;
>>F1=int（f,2,pi）
F1=
1/2*atan（1/2*2^（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.m
functiony=f（x）
y=1/sqrt（2*pi）*exp（-x.^2/2）;
保存后，在命令窗口键入
formatlong
>>[y,n]=quad（'f',-3,3）
则显示结果为：
y=
0.99729991863154
n=
57（表示被积函数f的调用次数）
3.2.3多重积分
3.2.3.1二重积分
一种方法：
由于积分区域的多样性和复杂性，通过人为地结合积分区域的图形，把二重积分转化为二次积分，用int命令进行计算；另一种方法：
利用Matlab提供的函数dblquad和quad2dggen计算，其调用格式如下：
·二重数值积分（在矩形区域上），使用dblquad函数，格式为：
y=dblquad（‘f’,x_m,x_M,y_m,y_M,tol）
·一般二重积分（在数值积分工具箱（NIT，即NumericalIntegrationToolbox）中提供了quad2dggen函数，该函数可从网上下载，网址：
fttp:
//
gration/nit.zip，安装之后，将其路径加载到Matlab路径下），格式为：
I=quad2dggen（被积函数名，下限函数名，上限函数名，y_m，y_M，tol）
例3-10计算
其中D为直线y=x和抛物线y=x2所围部分。
解：
由数学方法可得：
（1）画出积分区域示意图
>>symsxy
>>f=（2-x-y）/2;
>>y1=x;
>>y2=x^2;
>>ezplot（y1）;holdon
>>ezplot（y2）;
>>axis（[0,2,0,2]）
（2）确定积分限
a=fzero（'x-x^2',0）
a=
0
>>b=fzero（'x-x^2',1）
b=
1
这是第2次积分的上下限。
（3）积分运算
f_dy=int（f,y,x^2,x）%先对y积分
f_dy=
x-5/4*x^2-1/2*x*（x-x^2）+1/4*x^4
>>I=int（f_dy,a,b）%再对x积分
I=
11/120
例3-11试求下面的二次积分
解：
先定义函数，函数名为：
my2dfun.m
functionz=my2dfun（x,y）
globalkk;kk=kk+1;%定义全局变量，测定被积函数调用次数
z=exp（-x.^2/2）.*sin（x.^2+y）;
保存后在命令窗口执行下列命令：
clear
>>globalkk;kk=0;
>>y=dblquad（'my2dfun',-2,2,-1,1）,kk
y=
1.57449318974494
kk=
1038（表示被积函数调用次数）
3.2.3.2三重积分
三重积分的过程与二重积分过程相同，也是把三重积分转化成三次积分来计算，只不过在确定积分限时更繁琐而已。
在确定积分限时一般要结合三维的积分区域，所以一般先画出积分区域。
例3-12计算
其中区域V为三个坐标面及平面：
所围成的闭区域。
解：
由数学方法可得：
（1）画出积分区域示意图。
由于图形比较简单，所以我们没有画出。
（并且Matlab没有画三维符号函数的命令）
symsxyz
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/48
3.2.4Taylor展式
数学定理表明：
若f（x）在x=0点附近有直到n+1阶连续导数，那么：
其中：
实际上是要将函数f（x）表示成xn（n从0到无穷）的和的形式。
这完全可以用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=ex展开到6阶；
（2）把y=lnx在x=1点展开到6阶；
（3）把y=sinx在x=π/2点展开到6阶；
（4）把y=xt关于变量t展开到3阶。
解：
symsxt
y1=taylor（exp（-x））
y2=taylor（log（x）,6,1）
y3=taylor（sin（x）,pi/2,6）
y4=taylor（x^t,3,t）
运行结果为：
y1=
1-x+1/2*x^2-1/6*x^3+1/24*x^4-1/120*x^5
y2=
x-1-1/2*（x-1）^2+1/3*（x-1）^3-1/4*（x-1）^4+1/5*（x-1）^5
y3=
1-1/2*（x-1/2*pi）^2+1/24*（x-1/2*pi）^4
y4=
1+log（x）*t+1/2*log（x）^2*t^2
3.2.5Fourier级数
数学定理表明：
设函数f（x）已经展开为全区间上的一致收敛的三角级数：
其中：
虽然Matlab中还没有提供一个专门的命令用于求解Fourier级数，但根据上面的数学公式，也可以编写一个函数文件sfour.m，用于求函数f在区间
上的Fourier级数的系数：
function[a0,ak,bk]=sfour（f）
symsxn
a0=int（f,-pi,pi）/pi
an=int（f*cos（n*x）,-pi,pi）/pi;
an=simple（an）
bn=int（f*sin（n*x）,-pi,pi）/pi;
bn=simple（bn）
有了上面的函数，我们可以求得一些函数的Fourier级数的系数。
例3-14求函数y=x的Fourier级数的系数
解：
调用上面定义的函数sfour来计算：
>>symsx
>>f=x;
>>[a0,ak,bk]=sfour（f）
a0=
0
an=
0
bn=
2/n^2/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（xp,yp,zp）是等高面u=u（x,y,z）=C上的任一点，则称向量
为函数u=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方向上的数值梯度；
（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）%等高线图
>>holdon
>>quiver（px,py）%向量图
运行结果为：
数值梯度图
2.函数梯度和方向导数jacobian。
使用形式：
jacobian（f,v）
参数说明：
f：
函数向量或标量，当f为标量时，jacobian（f,v）=gradient（f）;
v：
自变量向量或者单个变量。
例3-16求：
（1）u=xyz在点播M（1,1,1）处的梯度；
（2）u=x2+2y2+3z2+xy+3x-2y-6z在点O（0,0,0）及A（1,2,3）处的梯度大小。
symsxyz
>>u1=x*y*z;
>>u2=x^2+2*y^2+3*z^2+x*y+3*x-2*y-6*z;
>>v=[x,y,z];
>>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=
111
>>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=
7712
3.3最小二乘法
在生产实践中，常常需要根据实际测量得到的一系列数据找出函数关系，通常叫做配曲线或找经验公式。
本节我们介绍一种找直线的方法，它是广泛使用的一种处理数据的方法。
例3-17在某实验中测得数据如下：
X
104
180
190
177
147
134
150
191
204
121
Y
100
200
210
185
155
135
170
205
235
125
由此要推出x和y的函数关系：
y=f（x）。
解：
先把这些数据点描绘出来，如图所示，观察x和y大概满足的函数关系。
先把x和y进行适当的调整，使自变量x的值从小到大排列。
>>x=[104180190177147134150191204121];
>>y=[100200210185155135170205235125];
>>[x,i]=sort（x）
x=
104121134147150177180190191204
i=
11065742389
>>y=y（i）
y=
100125135155170185200210205235
>>plot（x,y,'r*'）
>>holdon
这些数据点大致分布在一条直线上，即x和y有线性函数关系：
y=ax+b，由数学推断过程得：
，
其中n为
升级会员 