实验2 函数的极限、导数、积分与M文件的应用实验目的：
1、理解极限的概念；
2、掌握用matlab软件求函数极限的方法；
3、加深理解积分理论中分割、近似、求和、取极限的思想方法；
4、学习并掌握用matlab求导数，不定积分、定积分、二重积分、曲线积分的方法。
实验内容：
1、学习matlab命令：
应用matlab求极限的命令：
对应数学运算
matlab命令
limit(f)
limit(f,x,a)或者limit(f,a)
limit(f,x,a,’left’)
limit(f,x,a,’right’)
应用matlab求解方程的函数solve调用格式：
solve(f(x)),表示求方程f(x)=0的根。
理解极限概念,参照课本
M文件的创建与运行：
MATLAB是一个功能强大的程序设计语言，有着良好的交互计算的工作环境。包含MATLAB语言代码的文件称为M－文件，你可以使用文本编辑器来建立M－文件，然后可以象使用其它的MA TLAB函数或命令一样使用它们。
命令M文件及其运行：
命令文件没有输入参数，也不返回输出参数，只是一些命令行的组合。命令M文件中的命令可以访问MATLAB工作区中的所有变量，而且其中的所有变量也成为工作区的一部分，命令文件运行结束，命令文件产生的变量保留在工作区，直到关闭MATLAB或用命令删除。
（2）函数M文件及其调用：
在MATLAB编辑窗口还可建立函数M文件，我们可以根据需要建立自己的函数文件，它们能够像库函数一样方便地调用，从而可扩展MATLAB的功能。如果对于一类特殊的问题，建立起许多函数M文件，就能形成工具箱。格式如下：
function [输出参数列表]=函数名（形式参数列表）
函数体语句；
end
4、x=sym('x') 建立符号变量x
syms x y z 建立多个符号变量x,y,z
diff(f(x)) 求f(x)的一阶导数
diff(f(x),n) 求f(x)的n阶导数
diff(f(x,y),x) 求二元函数f(x,y)的关于x的偏导数
diff(f(x,y),x,n) 求二元函数f(x,y)的关于x的n阶偏导数
sum(x) x为向量时给出x的各个元素之和，若x为矩阵时给出x的列和
symsum(s,n) 表示求
symsum(s,k,m,n) 表示求
int(f(x)) 计算不定积分
int(f(x,y),x) 计算不定积分
int(f(x),a,b) 计算定积分
int(f(x,y),x,a,b) 计算定积分
观察数列当时的变化趋势。
解：输入matlab命令：
n=1:100;
xn=n./(n+1);
for i=1:100;
plot(n(i),xn(i),'b')
hold on
end
其中 for end是循环语句，循环体内的语句被执行100次，n(i)
表示n的第i个分量。由图可以看出，随着n的增大，点列与直线
无限接近，于是有结论：。
分析函数当时的变化趋势。
解：输入matlab命令
x=-1:0.01:1;
y=x.*sin(1./x);
plot(x,y)
hold on
plot(x,x,x,-x)
从图上看，f(x)随着|x|的减小，振幅越来越接近零，频率越来越高，做无限振荡。
分析函数当时的变化趋势。
解：输入matlab命令
x=-1:0.01:1;
y=sin(1./x);
plot(x,y)
从图上看，当时，在-1和1之间无限次振荡，极限不存在。仔细观察图形，回答为什么图像的有些峰值不是-1和1？
分析函数当时的变化趋势。
解：输入matlab命令
x=1:100;
y=(1+1./x).^x;
plot(x,y)
从图上可以看出当时，函数值与某常数无限接近，我们知道这个常数就是 e 。
例5 求极限。
解：输入matlab命令
syms x
limit((tan(x)-sin(x))/(sin(x)^3))
例6 解方程。
解：输入matlab命令
syms a b c x
f=a*x^2+b*x+c;
solve(f)
例7 解方程。
解：输入malab命令
syms x
f=x^5-5*x-1;
solve(f)
例8 分别建立命令式M文件和函数式M文件求向量的平均值和标准差。
解：命令式函数文件：
在文件编辑窗口输入如下命令：
x=rand(20);
n=length(x);
mean=sum(x)/n
stdev=sqrt(sum(x.^2)/n-mean.^2)
随便以一个文件名命名即可，例如保存为programme.m，在命令窗口输入programme即可得到结果。
函数式函数文件：
在文件编辑窗口输入如下命令：
function [mean,stdev]=stat(x)
n=length(x);
mean=sum(x)/n
stdev=sqrt(sum(x.^2)/n-mean.^2)
保存为文件stat.m，在命令窗口输入：
x=rand(20);
[mean,stdev]=stat(x)
即可得到结果。
例9 画出在（）处的切线及若干条割线，观察割线的变化趋势。
解：在曲线上另取一点，则的方程是：
，取，分别作出几条割线。
然后作出在处的切线。
h=[3,2,1,0.1,0.01];
a=(exp(h)-1)./h;
x=-1:0.1:3;
plot(x,exp(x),'r');
hold on
for i=1:5
plot(h(i),exp(h(i)),'r.')
plot(x,a(i)*x+1)
end
axis square
plot(x,x+1,'r')
从图上看出割线越来越接近切线。
例10 求隐函数所确定函数的导数。
解：，先求，再求,
则，于是输入matlab命令
syms x y
diffx=diff(log(x)+exp(-y/x)-exp(1),x);
diffy=diff(log(x)+exp(-y/x)-exp(1),y);
-diffx/diffy
例11 设隐函数，求。
解：利用matlab命令中的jacobian函数简化步骤，
syms x y z
a=jacobian(sin(x*y)+cos(y*z)+tan(z*x),[x,y,z]);
dz_dx=-a(1)/a(3)
dz_dy=-a(2)/a(3)
例12 利用定积分定义计算。
解：等距划分为20个子区间左矩形法求和近似：
x=linspace(0,1,21);
y=exp(x);
y1=y(1:20);
s1=sum(y1)/20
求得s1 =1.6757，图形演示如下
plot(x,y)
hold on
for i=1:20
fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i),y(i),0],'b')
end
右矩形法求和近似：
x=linspace(0,1,21);
y=exp(x);
y2=y(2:21);
s2=sum(y2)/20
求得s2 =1.7616，图形演示如下
plot(x,y)
hold on
for i=1:20
fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i+1),y(i+1),0],'b')
end
s1 和s2都可以作为该定积分的近似值。
直接用积分命令int也可以求，命令如下
syms x
int(exp(x),0,1)
例13 判断广义积分，,的敛散性，收敛时计算积分值。
解：对第一个积分：
syms p
int(1/x^p,x,1,inf)
由结果可以看出，当时，为无穷，当时，。
对第二个积分：
int(1/(2*pi)^(1/2)*exp(-x^2/2),-inf,inf)
对第三个积分：
int(1/(1-x)^2,0,2)
结果为，所以广义积分不收敛。
练习：
1、生物学中，有一个刻画生物群体的个体总量增长情况的著名方程——Logistic方程：Pn+1=kPn(1-Pn)，其中Pn为某一生物群体的第n代的个体总量与该群体所能达到的最大个体总量之比，,k为比例系数。
给定初值P0和比例系数k的值后，有方程就能生成一个数列，P0，P1，P2，……，Pn，……
生物学家为了预测群体总量的变化情况，就要研究这个数列，他们感兴趣的是：这个数列存在极限吗？数列中的项会出现周期性的变化情况吗？数列会不会出现无法预测的紊乱情况？（试给定不同的初值P0和比例系数k，研究生物群体变化情况）
2、完全数是指它的所有因子之和等于该完全数。6是一个完全数，因为6的因子为1，2，3且1+2+3=6,试分别建立命令式M文件和函数式M文件，找出10000以内的所有完全数，并对它做素因子分解。你能据此猜测完全数的通式吗？
3、绝对素数问题：一个素数，如果其中的任意两位数字交换后得到的数仍为素数，则将此素数成为绝对素数。例如：13和31都是素数，则13称为绝对素数。试利用matlab编程求出三位数中的所有绝对素数。要求分别使用命令式M文件和函数试M文件编写。
卫星轨道长度问题和覆盖面积问题：
人造地球卫星轨道可视为平面上的椭圆，我国第一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km，问该卫星的轨道长度？
一颗地球同步轨道通讯卫星的轨道位于地球的赤道平面内，且可近似认为是圆轨道，通讯卫星运行的角速度与地球自转的角速度相同，所以人们看到它在天空不动。问卫星距离地面的高度h应该为多少？并计算通讯卫星的覆盖面积。
5、路灯照明问题在一条20m宽的道路两侧，分别按装了一只2kW和3kW的路灯，它们离地面的高度分别为5m和6m。
在漆黑的夜晚，当两只路灯开启时，两只路灯连线的路面上最暗的点和最亮的点在那里？
如果3kW的路灯的高度可以在3m到9m之间变化，如何使路面上最暗点的亮度最大？作数值计算并绘图加以说明。
（注：，其中I是照度，r是光源到照射点的距离，是光源到照射点的光线与水平面的夹角，P是路灯的功率）
1、理解极限的概念；
2、掌握用matlab软件求函数极限的方法；
3、加深理解积分理论中分割、近似、求和、取极限的思想方法；
4、学习并掌握用matlab求导数，不定积分、定积分、二重积分、曲线积分的方法。
实验内容：
1、学习matlab命令：
应用matlab求极限的命令：
对应数学运算
matlab命令
limit(f)
limit(f,x,a)或者limit(f,a)
limit(f,x,a,’left’)
limit(f,x,a,’right’)
应用matlab求解方程的函数solve调用格式：
solve(f(x)),表示求方程f(x)=0的根。
理解极限概念,参照课本
M文件的创建与运行：
MATLAB是一个功能强大的程序设计语言，有着良好的交互计算的工作环境。包含MATLAB语言代码的文件称为M－文件，你可以使用文本编辑器来建立M－文件，然后可以象使用其它的MA TLAB函数或命令一样使用它们。
命令M文件及其运行：
命令文件没有输入参数，也不返回输出参数，只是一些命令行的组合。命令M文件中的命令可以访问MATLAB工作区中的所有变量，而且其中的所有变量也成为工作区的一部分，命令文件运行结束，命令文件产生的变量保留在工作区，直到关闭MATLAB或用命令删除。
（2）函数M文件及其调用：
在MATLAB编辑窗口还可建立函数M文件，我们可以根据需要建立自己的函数文件，它们能够像库函数一样方便地调用，从而可扩展MATLAB的功能。如果对于一类特殊的问题，建立起许多函数M文件，就能形成工具箱。格式如下：
function [输出参数列表]=函数名（形式参数列表）
函数体语句；
end
4、x=sym('x') 建立符号变量x
syms x y z 建立多个符号变量x,y,z
diff(f(x)) 求f(x)的一阶导数
diff(f(x),n) 求f(x)的n阶导数
diff(f(x,y),x) 求二元函数f(x,y)的关于x的偏导数
diff(f(x,y),x,n) 求二元函数f(x,y)的关于x的n阶偏导数
sum(x) x为向量时给出x的各个元素之和，若x为矩阵时给出x的列和
symsum(s,n) 表示求
symsum(s,k,m,n) 表示求
int(f(x)) 计算不定积分
int(f(x,y),x) 计算不定积分
int(f(x),a,b) 计算定积分
int(f(x,y),x,a,b) 计算定积分
观察数列当时的变化趋势。
解：输入matlab命令：
n=1:100;
xn=n./(n+1);
for i=1:100;
plot(n(i),xn(i),'b')
hold on
end
其中 for end是循环语句，循环体内的语句被执行100次，n(i)
表示n的第i个分量。由图可以看出，随着n的增大，点列与直线
无限接近，于是有结论：。
分析函数当时的变化趋势。
解：输入matlab命令
x=-1:0.01:1;
y=x.*sin(1./x);
plot(x,y)
hold on
plot(x,x,x,-x)
从图上看，f(x)随着|x|的减小，振幅越来越接近零，频率越来越高，做无限振荡。
分析函数当时的变化趋势。
解：输入matlab命令
x=-1:0.01:1;
y=sin(1./x);
plot(x,y)
从图上看，当时，在-1和1之间无限次振荡，极限不存在。仔细观察图形，回答为什么图像的有些峰值不是-1和1？
分析函数当时的变化趋势。
解：输入matlab命令
x=1:100;
y=(1+1./x).^x;
plot(x,y)
从图上可以看出当时，函数值与某常数无限接近，我们知道这个常数就是 e 。
例5 求极限。
解：输入matlab命令
syms x
limit((tan(x)-sin(x))/(sin(x)^3))
例6 解方程。
解：输入matlab命令
syms a b c x
f=a*x^2+b*x+c;
solve(f)
例7 解方程。
解：输入malab命令
syms x
f=x^5-5*x-1;
solve(f)
例8 分别建立命令式M文件和函数式M文件求向量的平均值和标准差。
解：命令式函数文件：
在文件编辑窗口输入如下命令：
x=rand(20);
n=length(x);
mean=sum(x)/n
stdev=sqrt(sum(x.^2)/n-mean.^2)
随便以一个文件名命名即可，例如保存为programme.m，在命令窗口输入programme即可得到结果。
函数式函数文件：
在文件编辑窗口输入如下命令：
function [mean,stdev]=stat(x)
n=length(x);
mean=sum(x)/n
stdev=sqrt(sum(x.^2)/n-mean.^2)
保存为文件stat.m，在命令窗口输入：
x=rand(20);
[mean,stdev]=stat(x)
即可得到结果。
例9 画出在（）处的切线及若干条割线，观察割线的变化趋势。
解：在曲线上另取一点，则的方程是：
，取，分别作出几条割线。
然后作出在处的切线。
h=[3,2,1,0.1,0.01];
a=(exp(h)-1)./h;
x=-1:0.1:3;
plot(x,exp(x),'r');
hold on
for i=1:5
plot(h(i),exp(h(i)),'r.')
plot(x,a(i)*x+1)
end
axis square
plot(x,x+1,'r')
从图上看出割线越来越接近切线。
例10 求隐函数所确定函数的导数。
解：，先求，再求,
则，于是输入matlab命令
syms x y
diffx=diff(log(x)+exp(-y/x)-exp(1),x);
diffy=diff(log(x)+exp(-y/x)-exp(1),y);
-diffx/diffy
例11 设隐函数，求。
解：利用matlab命令中的jacobian函数简化步骤，
syms x y z
a=jacobian(sin(x*y)+cos(y*z)+tan(z*x),[x,y,z]);
dz_dx=-a(1)/a(3)
dz_dy=-a(2)/a(3)
例12 利用定积分定义计算。
解：等距划分为20个子区间左矩形法求和近似：
x=linspace(0,1,21);
y=exp(x);
y1=y(1:20);
s1=sum(y1)/20
求得s1 =1.6757，图形演示如下
plot(x,y)
hold on
for i=1:20
fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i),y(i),0],'b')
end
右矩形法求和近似：
x=linspace(0,1,21);
y=exp(x);
y2=y(2:21);
s2=sum(y2)/20
求得s2 =1.7616，图形演示如下
plot(x,y)
hold on
for i=1:20
fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i+1),y(i+1),0],'b')
end
s1 和s2都可以作为该定积分的近似值。
直接用积分命令int也可以求，命令如下
syms x
int(exp(x),0,1)
例13 判断广义积分，,的敛散性，收敛时计算积分值。
解：对第一个积分：
syms p
int(1/x^p,x,1,inf)
由结果可以看出，当时，为无穷，当时，。
对第二个积分：
int(1/(2*pi)^(1/2)*exp(-x^2/2),-inf,inf)
对第三个积分：
int(1/(1-x)^2,0,2)
结果为，所以广义积分不收敛。
练习：
1、生物学中，有一个刻画生物群体的个体总量增长情况的著名方程——Logistic方程：Pn+1=kPn(1-Pn)，其中Pn为某一生物群体的第n代的个体总量与该群体所能达到的最大个体总量之比，,k为比例系数。
给定初值P0和比例系数k的值后，有方程就能生成一个数列，P0，P1，P2，……，Pn，……
生物学家为了预测群体总量的变化情况，就要研究这个数列，他们感兴趣的是：这个数列存在极限吗？数列中的项会出现周期性的变化情况吗？数列会不会出现无法预测的紊乱情况？（试给定不同的初值P0和比例系数k，研究生物群体变化情况）
2、完全数是指它的所有因子之和等于该完全数。6是一个完全数，因为6的因子为1，2，3且1+2+3=6,试分别建立命令式M文件和函数式M文件，找出10000以内的所有完全数，并对它做素因子分解。你能据此猜测完全数的通式吗？
3、绝对素数问题：一个素数，如果其中的任意两位数字交换后得到的数仍为素数，则将此素数成为绝对素数。例如：13和31都是素数，则13称为绝对素数。试利用matlab编程求出三位数中的所有绝对素数。要求分别使用命令式M文件和函数试M文件编写。
卫星轨道长度问题和覆盖面积问题：
人造地球卫星轨道可视为平面上的椭圆，我国第一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km，问该卫星的轨道长度？
一颗地球同步轨道通讯卫星的轨道位于地球的赤道平面内，且可近似认为是圆轨道，通讯卫星运行的角速度与地球自转的角速度相同，所以人们看到它在天空不动。问卫星距离地面的高度h应该为多少？并计算通讯卫星的覆盖面积。
5、路灯照明问题在一条20m宽的道路两侧，分别按装了一只2kW和3kW的路灯，它们离地面的高度分别为5m和6m。
在漆黑的夜晚，当两只路灯开启时，两只路灯连线的路面上最暗的点和最亮的点在那里？
如果3kW的路灯的高度可以在3m到9m之间变化，如何使路面上最暗点的亮度最大？作数值计算并绘图加以说明。
（注：，其中I是照度，r是光源到照射点的距离，是光源到照射点的光线与水平面的夹角，P是路灯的功率）
