数学实验课后习题解答.doc

数学实验课后习题解答
配套教材：王向东 戎海武 文翰 编著数学实验

王汝军编写

实验一　曲线绘图
【练习与思考】
画出下列常见曲线的图形。
以直角坐标方程表示的曲线：
1. 立方曲线
clear;
x=-2:0.1:2;
y=x.^3;
plot(x,y)

2. 立方抛物线
clear;
y=-2:0.1:2;
x=y.^3;
plot(x,y)
grid on

3. 高斯曲线
clear;
x=-3:0.1:3;

y=exp(-x.^2);
plot(x,y);
grid on
%axis equal

以参数方程表示的曲线
4. 奈尔抛物线
clear;
t=-3:0.05:3;
x=t.^3;y=t.^2;
plot(x,y)
axis equal
grid on

5. 半立方抛物线
clear;
t=-3:0.05:3;
x=t.^2;y=t.^3;
plot(x,y)
%axis equal
grid on

6. 迪卡尔曲线
clear;
a=3;t=-6:0.1:6;
x=3*a*t./(1t.^2);
y=3*a*t.^2./(1+t.^2);
plot(x,y)

7. 蔓叶线
clear;
a=3;t=-6:0.1:6;
x=3*a*t.^2./(1+t.^2);
y=3*a*t.^3./(1+t.^2);
plot(x,y)

8. 摆线
clear;clc;
a=1;b=1;
t=0:pi/50:6*pi;
x=a*(t-sin(t));
y=b*(1-cos(t));
plot(x,y);
axis equal
grid on

9. 内摆线（星形线）
clear;
a=1
t=0:pi/50:2*pi;
x=a*cos(t).^3;
y=a*sin(t).^3;
plot(x,y)

10. 圆的渐伸线（渐开线）
clear;
a=1;
t=0:pi/50:6*pi;
x=a*(cos(t)+t.*sin(t));
y=a*(sin(t)+t.*cos(t));
plot(x,y)
grid on

11. 空间螺线
clear
a=3;b=2;c=1;
t=0:pi/50:6*pi;
x=a*cos(t);
y=b*sin(t);
z=c*t;
plot3(x,y,z)
grid on

以极坐标方程表示的曲线：
12. 阿基米德线
clear;
a=1;
phy=0:pi/50:6*pi;
rho=a*phy;
polar(phy,rho,'r-*')

13. 对数螺线
clear;
a=0.1;
phy=0:pi/50:6*pi;
rho=exp(a*phy);
polar(phy,rho)

14. 双纽线
clear;
a=1;
phy=-pi/4:pi/50:pi/4;
rho=a*sqrt(cos(2*phy));
polar(phy,rho)
hold on

polar(phy,-rho)

15. 双纽线
clear;
a=1;
phy=0:pi/50:pi/2;
rho=a*sqrt(sin(2*phy));
polar(phy,rho)
hold on
polar(phy,-rho)

16. 四叶玫瑰线
clear;close
a=1;
phy=0:pi/50:2*pi;
rho=a*sin(2*phy);
polar(phy,rho)

17. 三叶玫瑰线
clear;close
a=1;
phy=0:pi/50:2*pi;
rho=a*sin(3*phy);
polar(phy,rho)

18. 三叶玫瑰线
clear;close
a=1;
phy=0:pi/50:2*pi;
rho=a*cos(3*phy);
polar(phy,rho)

实验二　极限与导数
【练习与思考】
1． 求下列各极限
（1） （2） （3）
clear;
syms n
y1=limit((1-1/n)^n,n,inf)
y2=limit((n^3+3^n)^(1/n),n,inf)
y3=limit(sqrt(n+2)-2*sqrt(n+1)+sqrt(n),n,inf)

y1 =1/exp(1)
y2 =3
y3 =0

（4） （5） （6）
clear;
syms x ;
y4=limit(2/(x^2-1)-1/(x-1),x,1)
y5=limit(x*cot(2*x),x,0)
y6=limit(sqrt(x^2+3*x)-x,x,inf)

y4 =-1/2
y5 =1/2
y6 =3/2

（7） （8） （9）
clear;
syms x m
y7=limit(cos(m/x),x,inf)
y8=limit(1/x-1/(exp(x)-1),x,1)
y9=limit(((1+x)^(1/3)-1)/x,x,0)

y7 =1
y8 =(exp(1) - 2)/(exp(1) - 1)
y9 =1/3

2． 考虑函数

作出图形，并说出大致单调区间；使用diff求，并求确切的单调区间。
clear;close;
syms x;
f=3*x^2*sin(x^3);
ezplot(f,[-2,2])
grid on
大致的单调增区间：[-2,-1.7],[-1.3,1.2],[1.7,2];
大致的单点减区间：[-1.7,-1.3],[1.2,1.7];

f1=diff(f,x,1)
ezplot(f1[-2,2])
line([-5,5],[0,0])
grid on
axis([-2.1,2.1,-60,120])

f1 =
6*x*sin(x^3) + 9*x^4*cos(x^3)

用fzero函数找的零点，即原函数的驻点
x1=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[-2,-1.7])
x2=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[-1.7,-1.5])
x3=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[-1.5,-1.1])
x4=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',0)
x5=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[1,1.5])
x6=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[1.5,1.7])
x7=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[1.7,2])

x1 =
-1.9948
x2 =
-1.6926
x3 =
-1.2401
x4 =
0
x5 =
1.2401
x6 =
1.6926
x7 =
1.9948
确切的单调增区间：[-1.9948,-1.6926],[-1.2401,1.2401],[1.6926,1.9948]
确切的单调减区间：[-2,-1.9948],[-1.6926,-1.2401],[1.2401,1.6926],[1.9948,2]
3． 对于下列函数完成下列工作，并写出总结报告，评论极值与导数的关系，
(i) 作出图形，观测所有的局部极大、局部极小和全局最大、全局最小值点的粗略位置；
(iI) 求所有零点（即的驻点）；
(iii) 求出驻点处的二阶导数值;
(iv) 用fmin求各极值点的确切位置;
(v) 局部极值点与有何关系?
(1)
(2)
(3)
clear;close;
syms x;
f=x^2*sin(x^2-x-2)
ezplot(f,[-2,2])
grid on

f =
x^2*sin(x^2 - x - 2)

局部极大值点为：-1.6,局部极小值点为为：-0.75，-1.6
全局最大值点为为：-1.6，全局最小值点为：-3
f1=diff(f,x,1)
ezplot(f1,[-2,2])
line([-5,5],[0,0])
grid on
axis([-2.1,2.1,-6,20])

f1 =
2*x*sin(x^2 - x - 2) + x^2*cos(x^2 - x - 2)*(2*x - 1)

用fzero函数找的零点，即原函数的驻点
x1=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[-2,-1.2])
x2=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[-1.2,-0.5])
x3=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[-0.5,1.2])
x4=fzero('2*x*sin(x^2-x2)+x^2*cos(x^2-x-2)*(2*x-1)',[1.2,2])

x1 =
-1.5326
x2 =
-0.7315
x3 =
-3.2754e-027
x4 =
1.5951
ff=@(x) x.^2.*sin(x.^2-x-2)
ff(-2),ff(x1),ff(x2),ff(x3),ff(x4),ff(2)

ff =
@(x)x.^2.*sin(x.^2-x-2)
ans =
-3.0272
ans =
2.2364
ans =
-0.3582
ans =
-9.7549e-054
ans =
-2.2080
ans =
0

实验三　级数
【练习与思考】
1. 用taylor命令观测函数的Maclaurin展开式的前几项, 然后在同一坐标系里作出函数和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向的图形的逼近的情况
(1)
clear;
syms x
y=asin(x);
y1=taylor(y,0,1)
y2=taylor(y,0,5)
y3=taylor(y,0,10)
y4=taylor(y,0,15)
x=-1:0.1:1;
y=subs(y,x);
y1=subs(y1,x);
y2=subs(y2,x);
y3=subs(y3,x);
y4=subs(y4,x);
plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3)

y1 =
0
y2 =
x^3/6 + x
y3 =
(35*x^9)/1152 + (5*x^7)/112 + (3*x^5)/40 + x^3/6 + x
y4 =
(231*x^13)/13312 + (63*x^11)/2816 + (35*x^9)/1152 + (5*x^7)/112 + (3*x^5)/40 + x^3/6 + x

(2)
clear;
syms x
y=atan(x);y1=taylor(y,0,3)
y2=taylor(y,0,5),y3=taylor(y,0,10),y4=taylor(y,0,15)
x=-1:0.1:1;
y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x);
y3=subs(y3,x);y4=subs(y4,x);
plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3)

y1 =
x
y2 =
x - x^3/3
y3 =
x^9/9 - x^7/7 + x^5/5 - x^3/3 + x
y4 =
x^13/13 - x^11/11 + x^9/9 - x^7/7 + x^5/5 - x^3/3 + x

(3)

clear;
syms x
y=exp(x^2);
y1=taylor(y,0,3)
y2=taylor(y,0,5)
y3=taylor(y,0,10)
y4=taylor(y,0,15)
x=-1:0.1:1;
y=subs(y,x);
y1=subs(y1,x);
y2=subs(y2,x);
y3=subs(y3,x);

y4=subs(y4,x);
plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3)

y1 =
x^2 + 1
y2 =
x^4/2 + x^2 + 1
y3 =
x^8/24 + x^6/6 + x^4/2 + x^2 + 1
y4 =
x^14/5040 + x^12/720 + x^10/120 + x^8/24 + x^6/6 + x^4/2 + x^2 + 1

(4)
clear;
syms x
y=sin(x)^2;
y1=taylor(y,0,1)
y2=taylor(y,0,5)
y3=taylor(y,0,10)
y4=taylor(y,0,15)
x=-pi:0.1:pi;
y=subs(y,x);
y1=subs(y1,x);
y2=subs(y2,x);
y3=subs(y3,x);
y4=subs(y4,x);
plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3)

y1 =
0
y2 =
x^2 - x^4/3
y3 =
- x^8/315 + (2*x^6)/45 - x^4/3 + x^2
y4 =
(4*x^14)/42567525 - (2*x^12)/467775 + (2*x^10)/14175 - x^8/315 + (2*x^6)/45 - x^4/3 + x^2

(5)
clear;
syms x
y=exp(x)/(1-x);
y1=taylor(y,0,3)
y2=taylor(y,0,5)
y3=taylor(y,0,10)
y4=taylor(y,0,15)
x=-1:0.1:0;
y=subs(y,x);
y1=subs(y1,x);
y2=subs(y2,x);
y3=subs(y3,x);
y4=subs(y4,x);
plotx,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3)

y1 =
(5*x^2)/2 + 2*x + 1
y2 =
(65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1
y3 =
(98641*x^9)/36288 + (109601*x^8)/40320 + (685*x^7)/252 + (1957*x^6)/720 + (163*x^5)/60 + (65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1
y4 =
(47395032961*x^14)/17435658240 + (8463398743*x^13)/3113510400 + (260412269*x^12)/95800320 + (13563139*x^11)/4989600 + (9864101*x^10)/3628800 + (98641*x^9)/36288 + (109601*x^8)/40320 + (685*x^7)/252 + (1957*x^6)/720 + (163*x^5)/60 + (65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1

(6)
clear;
syms x
ylog(x+sqrt(1+x^2));
y1=taylor(y,0,3)
y2=taylor(y,0,5)
y3=taylor(y,0,10)
y4=taylor(y,0,15)
x=-1:0.1:1;
y=subs(y,x);
y1=subs(y1,x);
y2=subs(y2,x);
y3=subs(y3,x);
y4=subs(y4,x);
plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3)

y1 =
x
y2 =
x - x^3/6
y3 =
(35*x^9)/1152 (5*x^7)/112 + (3*x^5)/40 - x^3/6 + x
y4 =
(231*x^13)/13312 - (63*x^11)/2816 + (35*x^9)/1152 - (5*x^7)/112 + (3*x^5)/40 - x^3/6 + x

2. 求公式中的数的值.
k=[4 5 6 7 8];
syms n
symsum(1./n.^(2*k),1,inf)

ans =
[ pi^8/9450, pi^10/93555, (691*pi^12)/638512875, (2*pi^14)/18243225, (3617*pi^16)/325641566250]

3. 利用公式来计算的近似值。精确到小数点后100位,这时应计算到这个无穷级数的前多少项?请说明你的理由.
解：Matlab代码为
clear;clc;close
epsl=1.0e-100;
ep=1;fn=1;a=1;n=1;
while ep>epsl
a=a+fn;
n=n+1;
fn=fn/n;
ep=fn;
end
fn
vpa(a,100)
n

fn =
8.3482e-101
ans =
2.71828182845904553488480814849026501178741455078125
n =
70
精确到小数点后100位,这时应计算到这个无穷级数的前71项,理由是误差小于10的负100次方，需要最后一项小于10的负100次方，由上述循环知n=70时最后一项小于10的负100次方，故应计算到这个无穷级数的前71项.
4. 用练习3中所用观测法判断下列级数的敛散性
(1)
clear;clc;
epsl=0.000001;
N=50000;p=1000;
syms n
Un=1/(n^2+n^3);
s1=symsum(Un,1,N);
s2=symsum(Un,1,N+p);
sa=vpa(s2-s1);
sa=setstr(sa)
sa=str2num(sa);
fprintf('级数')
disp(Un)
if sa
s1=symsum(Un,1,N);
s2=symsum(Un,1,N+p);
sa=vpa(s2-s1);
sa=setstr(sa);
sa=str2num(sa);
fprintf('级数')
disp(Un)
if saclear;clc;
epsl=0.00000000000001;
N=50000;p=100;
syms n
Un=1/sin(n);
s1=symsum(Un,1,N);
s2=symsum(Un,1,N+p);
sa=vpa(s2-s1);
sa=setstr(sa);
sa=str2num(sa);
fprintf('级数')
disp(Un)
if abs(sa)for k=1:100
s(k)=symsum(1/sin(n),1,k);
end
plot(s,'.')

发散
(4)
clear;clc;
epsl=0.0000001;
N=50000;p=1000;
syms n
Un=log(n)/(n^3);
s1=symsum(Un,1,N);
s2=symsum(Un,1,N+p);
sa=vpa(s2-s1);
sa=setstr(sa);
sa=str2num(sa);
fprintf('级数')
disp(Un)
if sa disp('发散')
end

级数log(n)/n^3收敛
clear;close
syms n
s=[];
for k=1:100
s(k)=symsum(log(n)/n^3,1,k);
end
plot(s,'.')

(5)
clear;close
syms n
s=[];he=0;
for k=1:100
he=he+factorial(k)/k^k;
s(k)=he;
end
plot(s,'.')

(6)
clear;clc;
epsl=0.0000001;
N=50000;p=1000;
syms n
Un=1/log(n)^n;
s1=symsum(Un,3,N);
s2=symsum(Un,3,N+p);
sa=vpa(s2-s1);
sa=setstr(sa);
sa=str2num(sa);
fprintf('级数')
disp(Un)
if sa
end
plot(s,'.')

(7)
clear;clc;
epsl=0.0000001;
N=50000;p=100;
syms n
Un=1/(log(n)*n);
s1=symsum(Un,3,N);
s2=symsum(Un,3,N+p);
sa=vpa(s2-s1);
sa=setstr(sa);
sa=str2num(sa);
fprintf('级数')
disp(Un)
if (sa)))发散

clear;close
syms n
s=[];
for k=3:300
s(k)=symsum(1/(n*log(n)),2,k);
end
plot(s,'.')

(8)
clear;clc;
epsl=0.0000001;
N=50000;p=100;
syms n
Un=(-1)^n*n/(n^2+1);
s1=symsum(Un,3,N);
s2=symsum(Un,3,N+p);
sa=vpa(s2-s1);
sa=setstr(sa);
sa=str2num(sa);
fprintf('级数')

disp(Un)
if (sa)
y2=1/(1+cos(x));
y3=1/(exp(x)+1);
y4=asin(x);
y5=sec(x)^3;
f1=int(y1)
f2=int(y2)
f3=int(y3)
f4=int(y4)
f5=int(y5)
dy=simplify(diff([f1;f2;f3;f4;f5]))

dy =
x*sin(x^2)
tan(x/2)^2/2 + 1/2
1/(exp(x) + 1)
asin(x)
(cot(pi/4 + x/2)*(tan(pi/4 + x/2)^2/2 + 1/2))/2 + 1/(2*cos(x)) + tan(x)^2/cos(x)

f1 =
-cos(x^2)/2
f2 =
tan(x/2)
f3 =
x - log(exp(x) + 1)
f4 =
x*asin(x) + (1 - x^2)^(1/2)
f5 =
log(tan(pi/4 + x/2))/2 + tan(x)/(2*cos(x))

2．(定积分)用trapz,quad,int计算下列定积分
,,,
解：Matlab代码为
clear;
x=(0+eps):0.05:1;
y1=sin(x)./x;
f1=trapz(x,y1)

f1 =0.9460
fun1=@(x)sin(x)./x;
f12=quad(fun1,0+eps,1)

f12 = 0.9461

f13=vpa(int('sin(x)/x',0,1),5)

f13 =0.94608

3．(椭圆的周长) 用定积分的方法计算椭圆的周长
解：椭圆的参数方程为
由参数曲线的弧长公式得

Matlab代码为
s=vpa(int('sqrt(5*sin(t)^2+4)','t',0,2*pi),5)

s =
15.865

4．（二重积分）计算数值积分
解：fxy=@(x,y)1+x+y;ylow=@(x)1-sqrt(1-x.^2);yup=@(x)1+sqrt(1-x.^2);
s=quad2d(fxy,-1,1,ylow,yup)

s =
6.2832
或符号积分法：
syms x y
xi=int(1+x+y,y,1-sqrt(1-x^2),1+sqrt(1-x^2));
s=int(xi,x,-1,1)
s =
2*pi
5．（假奇异积分）用trapz,quad8计算积分，会出现什么问题？分析原因，并求出正确的解。
解：Matlab代码为
clear
x=-1:0.05:1;
y=x.^(1/3).*cos(x);
s1=trapz(x,y)
fun5=@(x)x.^(1/3).*cos(x);

s2=quad(fun5,-1,1)
int('x^(1/3)*cos(x)','x',-1,1)

s1 =
0.9036 + 0.5217i
s2 =
0.9114 + 0.5262i
Warning: Explicit integral could not be found.
ans =
int(x^(1/3)*cos(x), x = -1..1) ，原函数不存在，不能用int函数运算。
用梯形法和辛普森法计算数值积分时，由于对负数的开三次方运算结果为复数，所以导致结果错误且为复数；
显然被积函数为奇函数，在对称区间上的积分等于0，此时可以这样处理：
（1）重新定义被积函数
%fun5.m
function y=fun5(x)
[m,n]=size(x);
for k=1:m
for l=1:n
y(k,l)=nthroot(x(k,l),3)*cos(x(k,l));
end
end
end
用辛普森法：
s=quad('fun5',-1,1)

s =
0

用梯形法
clear;
x=-1:0.01:1;
y=fun5(x);
s=trapz(x,y)

s =
-1.3878e-017

6．（假收敛现象）考虑积分，
（1）用解析法求;
clear;
syms x k;
Ik=int(abs(sin(x)),0,k*pi)

Warning: Explicit integral could not be found.
Ik =
int(abs(sin(x)), x = 0..pi*k)

（2）分别用trapz,quad和quad8求和,发现什么问题？
clear;
for k=4:2:8;
x=0:pi/1000:k*pi;
y=abs(sin(x));
trapz(x,y)
end

ans =
8.0000
ans =
12.0000
ans =
16.0000
for k=4:2:8
fun6=@(x)abs(sin(x));
quad(fun6,0,k*pi)
end

ans =
8.0000
ans =
12.0000
ans =
16.0000

7．(Simpson积分法)编制一个定步长Simpson法数值积分程序.计算公式为

其中为偶数,
解：Matlab代码为
%fun7.m
function y=fun7(f_name,a,b,n)
%f_name为被积函数
%[a,b]为积分区间
%n为偶数，用来确定步长h=(b-a)/n
if mod(n,2)~=0
disp('n必须为偶数')
return;
end
if nargin<4
n=100;
end
if nargin<3
disp('请输入积分区间')
end
if nargin==0
disp('error')
end
h=(b-a)/n;
x=a:h:b;
s=0;
for k=1:n+1
if k==1||k==(n+1)
xishu=1;
elseif mod(k,2)==0
xishu=4;
else
xishu=2;
end
s=s+feval(f_name,x(k))*xishu;
end
y=s*h/3;
end
8．(广义积分)计算广义积分
,,
并验证公式
.
解：Matlab代码为
clear;
syms x
s1=vpa(int(exp(-x^2)/(1+x^4),-inf,inf),5)
s2=quad(@(x)tan(x)./sqrt(x),0,1)
s3=quad(@(x)sin(x)./sqrt(1-x.^2),0,1)
s4=vpa(int(exp(-x^2/2)/sqrt(2*pi),-inf,inf),5)
s5=int(sin(x)./x,0+eps,inf)

s1 =
1.4348
s2 =
0.7968
s3 =
0.8933
s4 =
1.0
s5 =
pi/2 - sinint(1/4503599627370496)

实验五　二元函数的图形

【练习与思考】
1. 画出空间曲线在范围内的图形，并画出相应的等高线。
clear;
x=-30:0.5:30;y=-30:0.5:30;
[X,Y]=meshgrid(x,y);
Z=10*sin(sqrt(X.^2+Y.^2))./sqrt(1+X.^2+Y.^2);
mesh(X,Y,Z)

2. 根据给定的参数方程，绘制下列曲面的图形。
a) 椭球面
clear;
u=0:pi/50:2*pi;
v=0:pi/50:pi;
[U,V]=meshgrid(u,v);
x=3*cos(U).*sin(V);
y=2*cos(U).*cos(V);
z=sin(U);
mesh(x,y,z)

b) 椭圆抛物面
clear;
u=0:pi/50:pi/4;
v=0:pi/50:2*pi;
[U,V]=meshgrid(u,v);
x=3*U.*sin(V);
y=2*U.*cos(V);
z=4*U.^2;
mesh(x,y,z)
axis equal

c) 单叶双曲面
cle