程序为：
fplot（'sin（x）./x',[-1010-0.051]）,gtext（'sinx/x'）
图像为：
分析：
若改变其区间为fplot（'sin（x）./x',[-2030-0.43]）,gtext（'sinx/x'）
则图像为:
例3：
做函数x=cos3t,y=sin2t（0
t
2
）的图形;
程序为;
n=50
t=0:
pi/n:
2*pi;
x=cos（3*t）;y=sin（2*t）;
plot（x,y）
图像为：
结果分析：
以上两个例子表明自变量取值决定函数图像。
如n分别取1020304050等值时函数图像有很大变化。
例4：
在同一坐标系中作y=sin2xy=cosx（0
x2
）的图形
程序为：
x=0:
pi/15:
2*pi;
y1=sin（2*x）;
y2=cos（x）;
plot（x,y1,x,y2）
图像为:
例5
作分段函数
y=x+1,x<0;y=x,x
0的图形
程序;
x1=-5:
1/50:
0;
y1=x1+1;
x2=0:
1/50:
5;
y2=x2;
plot（x1,y1,x2,y2）,gtext（'y=x+1'）,gtext（'y=x'）
xlabel（'x'）
ylabel（'y'）
图像
例6
作函数ｒ＝１＋ｃｏｓｔ表示的心脏图形
程序
clear;
t=0:
2*pi/30:
2*pi;
r=1+cos（t）;
polar（t,r）
图形为：
例７
绘制极坐标下的曲线
ｐ＝acos（b+n
）的图形，并讨论a,b,n,取不同值的情形
程序为：
theta=0:
0.1:
2*pi;
fori=1:
4
a（i）=input（'a='）;
b（i）=input（'b='）;
n（i）=input（'n='）;
rho（i,:
）=a（i）*cos（b（i）+n（i）*theta）;
subplot（2,2,i）
polar（theta,rho（i,:
））;
end
当a=2，b=0,n=2,3,不同的取值运行后图形都不相同
分析：
以上为函数的平面图形的基本做法，若改变区间或变量将会得出不同的函数图像。
二.空间图形
例1作空间曲面
的图形
程序为：
x=-7.5:
0.5:
7.5;
y=x;
[X,Y]=meshgrid（x,y）;
R=sqrt（X.^4+Y.^4）+eps;
Z=sin（R）./R;
mesh（X,Y,Z）
xlabel（'x'）;
ylabel（'y'）;
zlabel（'z'）;
gridon
boxon
执行结果为：
例2
作函数（空间曲线）x=sint,y=cost,z=t,（0
t
15
）的图形
程序为：
t=-0:
pi/50:
15*pi;
x=sin（t）;
y=cos（t）;
z=t;
plot3（x,y,z）
xlabel（'x'）;ylabel（'y'）;zlabel（'z'）;
gridon
boxon
结果为：
例3.作椭球面
程序为：
ezmesh（'5*sin（u）*sin（v）','3*cos（u）*sin（v）','2*cos（v）',[0,2*pi,-pi*pi]）
结果为：
例4作双叶双曲面
程序为：
u=[0:
0.01:
0.4*pi0.6*pi:
0.01:
pi];
v=0:
0.01:
2*pi;
[u,v]=meshgrid（u,v）;
x=tan（u）.*cos（v）;
y=tan（u）.*sin（v）;
z=sec（u）;
mesh（x,y,z）
xlabel（'x'）;
ylabel（'y'）;
zlabel（'z'）;
图形为：
分析：
空间图形函数的做法和平面图形的做法有很多类似的地方，注意区间和变量的变化导致函数图像的变化。
二．函数的导数与最值
例1.铁路线上AB直线段长100Km，工厂C到铁路线上A处20Km，现在要在AB上选一点D，从D向C修一条直线公路，已知铁路运输每吨公里与公路运输每吨公里的运费之比为3:
5，为了是使原材料从B处运到工厂C的运费最省，D应在何处？
解：
设D点取在离A点距离x千米处,则BD=100-x,
CD=sqrt（AC^2+AD^2）=sqrt（400+x^2）
已知铁路运输每吨公里与公路运输每吨公里的运费之比为3:
5，所以设铁路上每吨千米运费为3元,公路每吨千米的运费为5元,设由B到C的运费为y元,故Y=3*（100-X）+5*sqrt（400+x^2）（0≤x≤100）
MATALAB程序为：
xmin=fminbnd（’3*（100-x）+5*sqrt（400+x^2）’,0,100）函数y在[0,100]上的极小点
xmin=15.0000
x=15;
y=3*（100-x）+5*sqrt（400+x^2）函数y在[0,100]上的极小值（运费）为y=380;
当x=0时，D点在A处的运费为y=400；
当x=100时，D点在B处的运费y=509.9020
综上所述,当x=15时,y取得最小值,即总运费最省.
因此,当D点取在与A相距15千米处时,总运费最省.
例2.求函数z=x^4+y^4-4xy+1的极值，并对图形进行观察。
MATLAB命令：
先画出该函数的图像：
>>clear;
>>x=0:
0.01:
1;
>>y=0:
0.01:
1;
>>[X,Y]=meshgrid（x,y）;
>>Z=X.^4+Y.^4-4.*X.*Y+1;
>>surf（X,Y,Z）
函数的曲面图:
>>clear;
>>symsxy;
>>Z=x^4+y^4-4*x*y+1;
>>Zx=diff（Z,x,1）;
>>Zy=diff（Z,y,1）;
>>simplify（Zx）
>>simplify（Zy）
结果为：
ans=
4*x^3-4*y
ans=
4*y^3-4*x
clear;
[x,y]=solve（'4*x^3-4*y,4*y^3-4*x','x','y'）
结果为：
x=
0
i
-i
-1
1
（1/2-1/2*i）*2^（1/2）
（-1/2+1/2*i）*2^（1/2）
（1/2+1/2*i）*2^（1/2）
（-1/2-1/2*i）*2^（1/2）
y=
0
-i
i
-1
1
-1/2*2^（1/2）-1/2*i*2^（1/2）
1/2*2^（1/2）+1/2*i*2^（1/2）
-1/2*2^（1/2）+1/2*i*2^（1/2）
1/2*2^（1/2）-1/2*i*2^（1/2）
>>clear;
>>symsxy;
>>Z=x^4+y^4-4*x*y+1;
>>A=diff（Z,x,2）
>>B=diff（diff（Z,x）,y）
>>C=diff（Z,y,2）
结果为:
A=
12*x^2
B=
-4
C=
12*y^2
因此当x=1,y=1时有极小值为1。
分析;函数的导数与最值在生活实际中应用很广泛，用MATLAB能准确快速的求解此类问题。
三．级数与函数逼近（拟合）
例1.求级数的和；
（1）
（2）
程序及远行结果为：
（1）n=1:
10;
S1=sum（1./（n.*（n+1）.*（n+2）））;
formatlong
S1symsn;
S2=symsum（1/（n*（n+1）*（n+2））,1,10）
vpa（S2）
结果：
S1=
0.24621212121212
S2=
65/264
ans=
.24621212121212121212121212121212
（2）
symsn;
symsum（1/n^2,1,inf）
ans=
1/6*pi^2
x=zeros（1,200）;y=zeros（1,200）;
fori=1:
200
y（i）=eval（symsum（1/n^2,1,i））;x（i）=i;
end
plot（x,y）;
xlabel（'n'）;
ylabel（'级数的部分和s（n）'）
x=zeros（1,200）;y=zeros（1,200）;
fori=1:
200
y（i）=eval（symsum（（-1）^n/n^2,1,i））;x（i）=i;
end
plot（x,y）;
xlabel（'n'）;
ylabel（‘级数的部分和s（n）‘）
结果：
（2）ans=
1/6*pi^2
分析：
通过图形课知道级数的部分和随着n值的变化而变化
例2
写出函数y=sinx的幂级数展开式，并用图形考察幂级数部分和的逼近情况
程序为：
x0=-2*pi:
0.01:
2*pi;
y0=sin（x0）;
symsx;
y=sin（x）;
plot（x0,y0,'r--'）,axis（[-2*pi,2*pi,-1.5,1.5]）;
holdon
p=taylor（y,x,2）,y1=subs（p,x,x0）;
line（x0,y1）
xlabel（'x轴'）;ylabel（'y轴'）;
gtext（'sin（x）'）;gtext（'2阶泰勒展式'）
图形：
分析：
由图知，k的取值控制展开的阶数；在0点附近逼近效果较好。
总结：
由比较判别法和计算结果知级数
发散，则
发散。
由根式判别法知该极限是收敛的。
由图可知，在[-2*pi,2*pi]区间上，用泰勒展开式，正弦函数y=sin（x）时，级数越高，越精确，即误差越小。
当n值越大时，级数值越接近精确
第二部分：
极限与微积分问题
一、实验目的
1.加深对极限概念的理解，学会用MATLAB计算极限。
2.加深对定积分概念的理解；学会用MATLAB计算定积分（数值近似计算和符号计算）。
3.了解常微分方程的基本概念；了解常微分方程的解析解；了解常微分方程的数值解
4.了解求解数值微分的基本方法；了解误差分析和步长优化。
5.会使用差分公式求解数值微分
二、实验内容
1.用数值计算和图形展示相结合研究数列和函数极限。
用符号演算和数值方法计算数列和函数极限。
2.用数值计算和图形展示结合研究函数的积分随分割细度的变化趋势；用数值方法和符号演算法计算定积分。
3.常微分方程模型的建立及求解。
4.中心差分公式，前向微分和后向微分公式。
三、相关知识
1.数列极限的定义及几何意义。
2.符号极限的MATLAB命令
3.定积分的定义；定积分的几何意义。
4.定积分存在的两个充分条件：
（1）f（x）在[a,b]连续一定可积；
（2）f（x）在[a,b]内有有限个第一类间断点，其余点连续，则f（x）在[a,b]可积。
5.常微分方程的解析解和数值解。
6.数值微分基本概念和意义；求符号倒数，数值微分的MATLAB命令
四．实验过程（操作过程，运行结果及结论）
一．数列的极限
求出下列极限的值
（1）lim（n^3+5^n）^（1/n）n趋于无穷
symsnan
an=limit（（n^3+5^n）^（1/n）,n,inf）
ans=
5
（2）.lim（sqrt（n+3）-3*sqrt（n+1）+sqrt（n））n趋于无穷
symsman
an=limit（（n+3）^（1/2）-3*（n+1）^（1/2）+n^（1/2）,n,inf）
ans=
-Inf
（3）.lim（cos（m/n））^nn趋于无穷
symsmnan
an=limit（cos（（m/n）^n）,n,inf）
ans=
1
（4）.lim（exp^（1/n））n趋于无穷
symsnan
an=limit（exp（1/n）,inf）
ans=
1
（5）.lim（（1/x）*sin（1/x）x趋于无穷
symsxan
an=limit（（1/x）*sin（1/x）,inf）
ans=
0
（6）.lim（1/x-1/（exp（x）-1）x趋于无穷
symsxan
an=limit（（1/x）-1/（exp（x）-1）,x,1）
ans=
（exp
（1）-2）/（exp
（1）-1）
（7）.lim（sin（a*x）/sin（b*x）x趋于0
Symsxan
An=limit（sin（a*x）/sin（b*x）,x,0）
ans=
a/b
（8）.lim（（1-cos（x））/x*sin（x）x趋于0
Symsxan
an=limit（（1-cos（x））/（x*sin（x））,x,0）
ans=
1/2
总结;掌握极限的基本求法，并能用MATLAB软件计算极限，给我们解决数学问题带来了极大的便利。
二．定积分的定义与计算
1.用定义计算定积分
例1
求积分
程序如下：
clearall;
f=inline（'（1+x^2）^-1'）;a=0;b=1;
n=30;
x=[];
x
（1）=a;
fork=1:
8
x（n+1）=b;s=0;
fori=1:
n-1
x（i+1）=（i+rand（））*（b-a）/n;
end
fori=1:
n
dxi=x（i+1）-x（i）;
c=x（i）+dxi*rand（）;
s=s+f（c）*dxi;
end
fprintf（'n=%g,s=%g\n',n,s）;
n=n*3;
end
结果为：
n=30,s=0.784935
n=90,s=0.78563
n=270,s=0.7854
n=810,s=0.785384
n=2430,s=0.785399
n=7290,s=0.785398
n=21870,s=0.785398
n=65610,s=0.785398
分析：
程序中分割小区间的个数n的初值取为30循序一次放大三倍，为了尽快获得结果。
若重复运行程序多次所得结果都不相同，当n=7290时根据结果可判断该定积分是存在的其值约为0.785398.
2.利用MATLAB模拟定积分定义域几何意义
例2设f（x）=sinx,求
并从图形观察随着分割点的增多，积分是否越来越接近定积分的值。
程序为：
functions=djfdf（f,a,b,n）
close;h=（b-a）/n;s=0;
fori=1:
n
x
（1）=a+（i-1）*h;
x
（2）=a+i*h;
x（3）=x
（2）;x（4）=x
（1）;t=（x（3）+x（4））/2;y（3）=feval（f,t）;
y（4）=y（3）;s=s+h*y（3）;
fill（x,y,[001]*i/n）;holdon;
end
fplot（f,[a,b]）;
holdoff
clearall;
clc
forn=1:
20;
f=inline（'sin（x）'）;
djfdf（f,0,2*pi,n）
pause（3）
end
结果：
分割点为20的图形
分割点为100的图形
分割点为600的图形
分析：
通过MATLAB计算定积分，从图形可以看出分割点越多越接近积分值
3.用定积分计算定积分值的简化
例3用上积分和与下积分和讨论函数f（x）=（1+x^2）^（-1）在区间【0,1】上的可积性
程序
cleart;
f=inline（'（x^2+1）^-1'）;a=0;b=1;
s0=1;
s1=0;
n=30;
t=[];
whileabs（s0-s1）>10^-4
t
（1）=a;
t（n+1）=b;
fori=1:
n-1
t（i+1）=（i+rand（））*（b-a）/n;
end
s0=0;s1==0;
fori=1:
n
s0=s0+f（t（i））*（t（i+1）-t（i））;
s1=s1+f（t（i+1））*（t（i+1）-t（i））;
end
n=n*3;
end
fprintf（'%s%s%g\n','1/（1+x^2）','在[0，1]上积分的近似值为'，s0）；
结果
1/（1+x^2）在[0，1]上积分的近似值为0.785438
结果分析：
当定积分存在时，如当被积函数f（x）连续时，由于无论
怎样划分区间，也无论在每个小区间内怎样取点
i,只要
趋近于0，积分和都收敛于相同的极限，因此，为简化计算，通常采取n等分区间，幷统一取每个小区间的左端点或右端点为
i，这样，积分和就变为h
或h
其中h=（b-a）/n.
对比6.7和6.8知抛物线法比梯形法收敛要快
4.定积分的近似计算方法
例4
用梯形公式计算定积分4/（1+x^2）在[0,1]的积分值为pi.
程序
formatlong;
n=100;
a=0;
b=1;
intsum=0;
symsxfx
fx=4/（1+x^2）;
fori=1:
n
xj=a+（i-1）*（b-a）/n;
xi=a+i*（b-a）/n;
fxj=subs（fx,'x',xj）;
fxi=subs（fx,'x',xi）;
intsum=intsum+（fxj+fxi）*（b-a）/（2*n）;
end
intsum
integrate=int（fx,0,1）
integrate=double（integrate）
fprintf（'相对误差是：
%d\n\n',abs（（intsum-integrate）/integrate））
结果
intsum=
3.14157598692313
integrate=
pi
integrate=
3.14159265358979
相对误差是：
5.305165e-006
5.用MATLAB软件中的积分函数
例5求
程序为
>>x=0:
0.01:
1;
>>y=x.*exp（x）;
>>trapz（x,y）
结果：
ans=
1.00003697125446
6.定积分的应用举例
例6某机器使用t周后的转售价格函数为R（t）=0.75Aexp（-t/96）（元），其中A为机器最初的价格。
在任何时间t,机器开动就能产生的利润为L=0.25Aexp（-t/48）（元），试问机器使用了多长时间后转售出去，才能使总利润最大？
利润是多少？
机器卖了多少钱？
clearall
clc
symsxtA
R=3*A/4*exp（-t/96）;
L=A/4*exp（-t/48）
f=subs（R,t,x）+int（L,0,x）;
x0=solve（diff（f,x））
fx02=subs（diff（f,x,2）,x,x0）;
fx0=subs（f,x,x0）
Lmax=fx0-A
Rx0=subs（R,t,x0）
Vpa（x0,1）
结果L=
1/4*A*exp（-1/48*t）
x0=
96*log（32）
fx0=
3075/256*A
Lmax=
2819/256*A
Rx0=
3/128*A
ans=
333
三．常微分方程和人口模型
例：
计算1952——2006年55年玉溪的人口六年增长率
>>x=[9541951000414106423612935881508149...
16332161766938188562320167802106416];
>>fori=1:
5
r（i）=（x（i+1）-x（i））/x（i）;
end
>>r
>>x=[9541951000414106423612935881508149...
16332161766938188562320167802106416];
>>fori=1:
9
r（i）=（x（i+1）-x（i））/x（i）;
end
>>r
Logistic模型
>>x=[9541951000414106423612935881508149...
16332161766938188562320167802106416];
>>fori=1:
9
r（i）=（x（i+1）-x（i））/x（i）;
end
>>x1=x（2:
10）;
>>plot（x1,r）
模型的参数估计
functionf=xt（g）
t=0:
9;
x=[954195100041410642361293588150814916332161766938188562320167802106416];
f=x-g
（1）./（1+（g
（1）./954195-1）.*exp（-g
（2）.*t））
g0=[3300,0.2];
g=lsqnonlin（'xt
升级会员 