王向东数学实验课后习题解答(第二篇2.1-2.10)

数学实验课后习题解答
配套教材：王向东戎海武文翰编著数学实验
王汝军编写
实验一曲线绘图【练习与思考】
画出下列常见曲线的图形。

以直角坐标方程表示的曲线：
1.立方曲线
3
x y=
clear;
x=-2:0.1:2; y=x.^3; plot(x,y)
2.立方抛物线
3x y=
clear;
y=-2:0.1:2; x=y.^3; plot(x,y) grid on
3.高斯曲线
2
x
e y-=
clear;
x=-3:0.1:3;
y=exp(-x.^2); plot(x,y); grid on
%axis equal
以参数方程表示的曲线
4. 奈尔抛物线)(,3
223x y t y t x ===
clear;
t=-3:0.05:3; x=t.^3;y=t.^2; plot(x,y) axis equal grid on
5. 半立方抛物线2323,()x t y t y x ===
clear;
t=-3:0.05:3; x=t.^2;y=t.^3; plot(x,y) %axis equal grid on
6.
迪卡尔曲线233
22
33,(30)11at at x y x y axy t t
==+-=++ clear;
a=3;t=-6:0.1:6; x=3*a*t./(1+t.^2); y=3*a*t.^2./(1+t.^2); plot(x,y)
7.
蔓叶线2332
22
,()11at at x x y y t t a x
===++- 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. 摆线)cos 1(),sin (t b y t t a x -=-=
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. 内摆线（星形线）)(sin ,cos 3
2323233a y x t a y t a x =+==
clear;
a=1;
t=0:pi/50:2*pi; x=a*cos(t).^3; y=a*sin(t).^3; plot(x,y)
10. 圆的渐伸线（渐开线）)cos (sin ),sin (cos t t t a y t t t a x -=+=
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. 空间螺线ct z t b y t a x ===,sin ,cos
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. 阿基米德线0,≥=r a r
ϕ
clear; a=1;
phy=0:pi/50:6*pi; rho=a*phy;
polar(phy,rho,'r-*')
13. 对数螺线ϕ
a e r =
clear; a=0.1;
phy=0:pi/50:6*pi; rho=exp(a*phy); polar(phy,rho) 14. 双纽线))()((2cos 2222222
2
y x a y x a r -=+=ϕ
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. 双纽线)2)((2sin 22222
2
xy a y x a r =+=ϕ
clear; a=1;
phy=0:pi/50:pi/2;
rho=a*sqrt(sin(2*phy)); polar(phy,rho) hold on
polar(phy,-rho)
16. 四叶玫瑰线0,2sin ≥=r a r ϕ
clear;close a=1;
phy=0:pi/50:2*pi; rho=a*sin(2*phy); polar(phy,rho)
17. 三叶玫瑰线0,3sin ≥=r a r ϕ
clear;close a=1;
phy=0:pi/50:2*pi; rho=a*sin(3*phy); polar(phy,rho)
18. 三叶玫瑰线0,3cos ≥=r a r ϕ
clear;close a=1;
phy=0:pi/50:2*pi; rho=a*cos(3*phy); polar(phy,rho)
实验二 极限与导数
【练习与思考】
1． 求下列各极限
（1）n
n n
)11(lim -
∞
→ （2）n n
n n 3lim 3+∞→ （3）)122(
lim n n n n ++-+∞
→
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）)1
112(
lim 21
---→x x x （5）x x x 2cot lim 0→ （6）)3(lim 2
x x x x -+∞→
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）x x x m )(cos lim ∞→ （8）)1
1
1(lim 1--→x x e x （9）x x x 11lim
3
0-+→ 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． 考虑函数
22),sin(3)(32<<-=x x x x f
作出图形，并说出大致单调区间；使用diff 求)('x f ，并求)(x f 确切的单调区间。

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 函数找)('x f 的零点，即原函数)(x f 的驻点 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) 求
)('x f 所有零点（即)(x f 的驻点）； (iii) 求出驻点处)(x f 的二阶导数值;
(iv) 用fmin 求各极值点的确切位置; (v) 局部极值点与)("),('x f x f 有何关系?
(1) ]2,2[),2sin()(22-∈--=x x x x x f (2) ]3,3[,10203)(35-∈+-=x x x x f (3)
]3,0[,2)(23∈---=x x x x x f
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 =
用fzero 函数找)('x f 的零点，即原函数)(x f 的驻点
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-x-2)+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 命令观测函数
)(x f y =的Maclaurin 展开式的前几项, 然后在同一坐标系里作出函数
)(x f y =和它的Taylor 展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向)(x f y =的图形的逼近的情况
(1) x x f arcsin )(=
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
(2) x x f arctan )(= 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 =
(3)
2
)
(x e x
f
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) x x f 2
sin )(= 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
(5)
x
e x
f x
-=
1)(
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);
plot(x,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 +
(6) )1ln()(2
x x x f ++= clear; syms x
y=log(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. 求公式,,2,1,1212 ==∑∞
=k m n k
k
n k π中的数8,7,6,5,4,=k m k 的值. 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. 利用公式
e n n =∑∞
=0
!1
来计算e 的近似值。

精确到小数点后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)
∑∞
=+1321n n
n
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<epsl
disp('收敛') else
disp('发散') end
级数1/(n^3 + n^2)收敛 clear;close syms n s=[];
for k=1:100
s(k)=symsum(1/(n^3 + n^2),1,k); end
plot(s,'.')
(2)
∑∞
=12
1n n n
clear;clc;
epsl=0.000001; N=50000;p=1000; syms n
Un=1/(n*2^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 sa<epsl
disp('收敛') else
disp('发散') end
级数1/(2^n*n)收敛 clear;close syms n s=[];
for k=1:100
s(k)=symsum(1/(2^n*n),1,k); end
plot(s,'.')
(3)
∑∞
=1
1sin
n n
clear;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)<epsl disp('收敛') else
disp('发散') end
级数1/sin(n)发散
clear;close syms n s=[];
for k=1:100
s(k)=symsum(1/sin(n),1,k); end
plot(s,'.')
发散 (4)
3
1ln n n
n
∞
=∑
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<epsl
disp('收敛') else
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)
∑∞
=1!n n
n
n clear;close syms n s=[];he=0; for k=1:100
he=he+factorial(k)/k^k; s(k)=he; end
plot(s,'.')
(6)
∑∞
=3)
(ln 1n n n
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<epsl
disp('收敛') else
disp('发散') end
级数1/log(n)^n 收敛
clear;close syms n s=[];
for k=3:100
s(k)=symsum(1/log(n)^n,3,k); end
plot(s,'.')
(7)
∑∞
=1
ln 1
n n n 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)<epsl disp('收敛') else
disp('发散') end
级数1/(n*log(n)) 发散
clear;close syms n s=[];
for k=3:300
s(k)=symsum(1/(n*log(n)),2,k); end
plot(s,'.')
(8) ∑∞
=+-1
21)1(n n n n
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)<epsl disp('收敛') else
disp('发散')
end
级数((-1)^n*n)/(n^2 + 1)收敛
clear;close syms n s=[];
for k=3:300
s(k)=symsum((-1)^n*n/(n^2+1),2,k); end
plot(s,'.')
实验四 积分
【练习与思考】
1．(不定积分)用int 计算下列不定积分，并用diff 验证
⎰
dx x x 2
sin ，⎰+x
dx
cos 1，
⎰+1x e dx ，⎰
xdx arcsin ，⎰xdx 3
sec 解：Matlab 代码为：
syms x
y1=x*sin(x^2); 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 计算下列定积分
⎰1
0sin dx x x ,
⎰
1
dx x x
,⎰π
20
)2sin(dx x e x
,⎰-1
2
dx e x
解：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．(椭圆的周长) 用定积分的方法计算椭圆14
92
2=+y x 的周长
解：椭圆的参数方程为⎩
⎨⎧==t y t
x sin 2cos 3
由参数曲线的弧长公式得
22
s π
π
==⎰
⎰
20
π
=⎰
Matlab 代码为
s=vpa(int('sqrt(5*sin(t)^2+4)','t',0,2*pi),5)
s =
15.865
4．（二重积分）计算数值积分
⎰⎰≤+++y
y x dxdy y x 222
)1(
解：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计算积分⎰-1
13/1cos xdx x ，会出现什么问题？分析原因，并求出正确的解。

解：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 =
用梯形法 clear;
x=-1:0.01:1; y=fun5(x); s=trapz(x,y)
s =
-1.3878e-017
6．（假收敛现象）考虑积分⎰=π
k dx x k I 0
sin )
(，
（1）用解析法求)(k I ; 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求)6(),4(I I 和)8(I ,发现什么问题？ 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 法数值积分程序.计算公式为
)42424(3114321+-+++++++=
≈n n n n f f f f f f f h
S I 其中n 为偶数,.1,,2,1),)1((,+=-+=-=n i h i a f f n
a
b h i
解：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．(广义积分)计算广义积分
⎰∞
∞-+-dx x x 421)
exp(,⎰10)tan(dx x x ,⎰-1021sin dx x x
并验证公式
⎰∞
∞-=-,
12)2exp(2
dx x π
⎰
∞
=0
2
sin π
dx x x .
解：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.
画出空间曲线z
=
在30,30x y -<<范围内的图形，并画出相应的等高线。

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) 椭球面3cos sin ,
2cos cos ,sin x u v y u v z u ===
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) 椭圆抛物面23sin ,2cos ,4x u v y u v z u ===
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) 单叶双曲面3sec sin ,2sec cos ,4tan x u v y u v z u ===
clear;
u=0:pi/15:pi; v=0:pi/15:2*pi;
[U,V]=meshgrid(u,v); x=3*sec(U).*sin(V); y=2*sec(U).*cos(V); z=4*tan(U); mesh(x,y,z)
d)
双曲抛物面22
,,3
u v x u y v z -===
clear
u=-3:0.1:3;
[U,V]=meshgrid(u); x=U; y=V;
z=(U.^2-V.^2)/3; mesh(x,y,z)
e) 旋转面ln sin ,ln cos ,x u v y u v z u ===
clear; u=1:0.1:5;
v=0:pi/30:2*pi;
[U,V]=meshgrid(u,v); x=log(U).*sin(V); y=log(U).*cos(V); z=U; mesh(x,y,z) axis equal
f) 圆锥面sin ,
cos ,x u v y u v z u ===
clear;
u=-5:0.1:5; v=0:pi/30:2*pi;
[U,V]=meshgrid(u,v); x=(U).*sin(V); y=(U).*cos(V); z=U;
mesh(x,y,z) axis equal
g) 环面(30.4cos )cos ,(30.4cos )sin ,0.4sin x u v y u v z v =+=+=
clear;
u=0:pi/30:2*pi; v=u;
[U,V]=meshgrid(u,v); x=(3+0.4*cos(U)).*cos(V); y=(3+0.4*cos(U)).*sin(V); z=0.4*sin(V); mesh(x,y,z)
h) 正螺面sin ,
cos ,4x u v y u v z v ===
clear;
u=0:pi/30:pi; v=0:pi/30:10*pi;
[U,V]=meshgrid(u,v); x=U.*sin(V); y=U.*cos(V); z=4*V; mesh(x,y,z) colorbar
3. 在一丘陵地带测量搞程,x 和y 方向每隔100米测一个点,得高程见表5-2,试拟合一曲面,确定合适的模型,并由此找出最高点和该点的高程.
clc;clear;
x1=[100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400]; x2=[100 200 300 400 100 200 300 400 100 200 300 400 100 200 300 400]; y=[636 698 680 662 697 712 674 626 624 630 598 552 478 478 412 334]'; x=[x1',x2']; x0=[1 1 1 1 1];
beta=lsqcurvefit('heigh',x0,x,y) %绘图：
a1=100:5:400; a2=a1;
[xx1,xx2]=meshgrid(a1,a2);
Z=beta(1)+beta(2)*xx1+beta(3)*xx2+beta(4)*xx1.^2+beta(5)*xx2.^2; mesh(xx1,xx2,Z)
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the default value of the function tolerance.
beta =
Columns 1 through 3
538.4375 1.4901 0.6189
Columns 4 through 5
-0.0046 -0.0017
contour(xx1,xx2,Z,30),colorbar
%计算最高点及高程
x0=[100,100];
options=optimset('largescale','off');
%设置下界
lb=[0,0];
%无上界
ub=[];
[x,fval]=fmincon('height',x0,[],[],[],[],lb,ub,[],options)
Warning: Options LargeScale = 'off' and Algorithm = 'trust-region-
reflective' conflict.
Ignoring Algorithm and running active-set algorithm. To run trust-region-reflective, set
LargeScale = 'on'. To run active-set without this warning, use Algorithm = 'active-set'.
> In fmincon at 445
Local minimum possible. Constraints satisfied.
fmincon stopped because the predicted change in the objective function
is less than the default value of the function tolerance and constraints were satisfied to within the default value of the constraint tolerance.
No active inequalities.
x =
161.9676 182.0320
fval =
400400 450 500 550 600 650 700
heigh 和height 两个函数分别定义如下：(应写在m 文件中) %heigh.m
function f=heigh(beta,xdata) xx1=xdata(:,1); xx2=xdata(:,2);
f=beta(1)+beta(2)*xx1+beta(3)*xx2+beta(4)*xx1.^2+beta(5)*xx2.^2; end
%height.m
function y=height(x)
y=-(538.4375+1.4901*x(1)+0.6189*x(2)-0.0046*x(1).^2-0.0017*x(2).^2); end
实验六 多元函数的极值
【练习与思考】 1.求1444+-+=xy y x z
的极值，并对图形进行观测。

解：Maltab 代码为 syms x y;
z=x^4+y^4-4*x*y+1; dzx=diff(z,x); dzy=diff(z,y);
s=solve(dzx,dzy,x,y); x=s.x.' y=s.y.'
x =
[ 0, 1, -1, (-1)^(3/4), -(-1)^(3/4), -i, i, -(-1)^(3/4)*i, (-1)^(3/4)*i] y =
[ 0, 1, -1, (-1)^(1/4), -(-1)^(1/4), i, -i, (-1)^(1/4)*i, -(-1)^(1/4)*i]
经计算可知，函数的驻点为(0,0)、(1,1)、(-1,-1) ezmeshc(z,[-2,2,-2,2])
从图形上观测可知，(1,1)、(-1,-1)为极值点，(0,0)不是极值点。

clear syms x y;
z=x^4+y^4-4*x*y+1; dzx=diff(z,x); A=diff(z,x,2) B=diff(dzx,y) C=diff(z,y,2)
A = 12*x^2
B = -4
C =
12*y^2
由判别法可知(1,1)、(-1,-1)均为极小值点。

2.求函数
()222,y x y x f +=在圆周122=+y x 的最大值和最小值。

解：构造Lagrange 函数
)1(2),(2222-+++=y x y x y x L λ
求Lagrange 函数的自由极值.先求L 关于λ,,y x 的一阶偏导数，再解正规方程可得所求的极值点，Matlab 代
码为 clear; syms x y k
L=x^2+2*y^2+k*(x^2+y^2-1); dlx=diff(L,x); dly=diff(L,y); dlk=diff(L,k);
s=solve(dlx,dly,dlk,x,y,k); k=s.k' x=s.x' y=s.y' k =
[ -1, -2, -1, -2] x =
[ 1, 0, -1, 0] y =
[ 0, 1, 0, -1] t=0:pi/50:2*pi; x=cos(t); y=sin(t);
z=x.^2+2*y.^2; plot3(x,y,z)
可得点（1,0）、（0,1）（-1,0）、（0，-1）为函数的条件极值点，经判断函数()222,y x y x f +=在
（1,0）、（-1,0）取得极小值，在（0,1）、（0，-1）取得极大值。

3.在球面1222
=++z y x 求出与点(3,1,-1)距离最近和最远点。

解：
设
球
面
上
的
点
为
(x,y,z),
则
此
点
与
点
(3,1,-1)
的
距
离
为
()222)1()1()3(,,++-+-=z y x z y x d 且(x,y,z)满足1222=++z y x ；构造
Lagrange 函
数
)1()1()1()3(),,,(222222-+++++-+-=z y x z y x z y x L λλ
求Lagrange 函数的自由极值.先求L 关于λ,,,z y x 的一阶偏导数，再解正规方程可得所求的极值点，Matlab
代码为
clear clear;
syms x y z k
L=(x-3)^2+(y-1)^2+(z+1)^2+k*(x^2+y^2+z^2-1);
dlx=diff(L,x); dly=diff(L,y); dlz=diff(L,z); dlk=diff(L,k);
s=solve(dlx,dly,dlz,dlk,x,y,z,k); x=s.x' y=s.y' z=s.z' k=s.k'
x =
[ (3*11^(1/2))/11, -(3*11^(1/2))/11] y =
[ 11^(1/2)/11, -11^(1/2)/11] z =
[ -11^(1/2)/11, 11^(1/2)/11] k =
[ 11^(1/2) - 1, - 11^(1/2) - 1]
vpa(eval(L),5)
ans =
[ 5.3668, 18.633]
得到条件极值点为(
,111111-、(111111
--，经判断，球面
1222=++z y x 上与点(3,1,-1)距离最近的点为(
)111111
-，最远的点
(111111
--。

4.求函数
z y x z y x f 32),,(++=在平面1=+-z y x 与柱面122=+y x 的交线上的最大值。

解：构造Lagrange 函数
)1()1(32),,,(2221-++-+-+++=y x z y x z y x z y x L λλλ
求Lagrange 函数的自由极值.先求L 关于21,,,,λλz y x 的一阶偏导数，再解正规方程可得所求的极值点，
Matlab 代码为 clear;
syms x y z k1 k2
L=x+2*y+3*z+k1*(x-y+z-1)+k2*(x^2+y^2-1); dlx=diff(L,x); dly=diff(L,y); dlz=diff(L,z); dlk1=diff(L,k1); dlk2=diff(L,k2);
s=solve(dlx,dly,dlz,dlk1,dlk2,x,y,z,k1,k2); x=s.x' y=s.y' z=s.z' k1=s.k1' k2=s.k2'
x =
[ (2*29^(1/2))/29, -(2*29^(1/2))/29] y =
[ -(5*29^(1/2))/29, (5*29^(1/2))/29] z =
[ 1 - (7*29^(1/2))/29, (7*29^(1/2))/29 + 1] k1 =
[ -3, -3] k2 =
[ 29^(1/2)/2, -29^(1/2)/2] eval(L)
ans =
[ 3 - 29^(1/2), 29^(1/2) + 3] 经判断可知，函数z y x z y x f 32),,(++=在平面1=+-z y x 与柱面122=+y x 的交线上的最
大值为3+
5.求函数2
2y x z +=在三条直线1,1,1=+==y x y x 所围区域上的最大值和最小值。

解：显然此函数的驻点为(0,0)不在此区域内，因此该函数的最大值和最小值点应在三条边界上，下面分别求此函数在这三条边界上的最大值和最小值，Matlab 代码如下 （1）求函数在直线边界x=1,0≤y ≤1上的最大值和最小值
将x=1代入原函数，则二元函数变为一元函数
z=1+y 2,0≤y ≤1
最大值点为y=1，最大值为2，最小值点为y=0，最小值为1； （2）求函数在直线边界y=1,0≤x ≤1上的最大值和最小值
将x=1代入原函数，则二元函数变为一元函数
z=1+x 2,0≤x ≤1
最大值点为x=1，最大值为2，最小值点为x=0，最小值为1； (3)求函数在直线边界x+y=1,0≤x ≤1上的最大值和最小值
将y=1-x 代入原函数，则二元函数变为一元函数
z=(1-x)2+x 2,0≤x ≤1
用Matlab 命令求此函数的最大和最小值点 先求驻点 clear; syms x
z=(1-x)^2+x^2; dzx=diff(z,x); x=solve(dzx,x) x =1/2
z1=eval(z)
计算在驻点处的函数值 z1 =1/2
计算在区间端点处的函数值 z2=subs(z,0) z3=subs(z,1) z2 = 1 z3 = 1
比较函数在各点处的函数值可知函数的最大值点为(1,1)，对应的最大值为2， 最小值点为(1/2,1/2)，最小值为1/2。

实验七 常微分方程
【练习与思考】
1．求下列微分方程的解析解
a) 一阶线性方程2'3
=-y x y dsolve('Dy-x^3*y=2','x')
ans =
C2*exp(x^4/4) + exp(x^4/4)*int(2/exp(x^4/4), x)
b) 贝努利方程
0'2=--y xy y
dsolve('Dy-x*y^2-y=0','x')
ans =
0 exp(x)/(C4 - exp(x)*(x - 1))
c) 高阶线性齐次方程02'3"'"=+--y
y y y
dsolve('D3y-D2y-3*Dy+2*y=0')
ans =
C8*exp(2*t) + C6*exp(t*(5^(1/2)/2 - 1/2)) + C7/exp(t*(5^(1/2)/2 + 1/2))
d) 高阶线性非齐次方程x y y y sin 32'3"=+- dsolve('D2y-3*Dy+2*y=3*sin(x)','x')
ans =
(9*cos(x))/10 + (3*sin(x))/10 + C11*exp(x) + C10*exp(2*x)
e) 欧拉方程223
4'3"'"x xy y x y x
=-+
dsolve('x^3*D3y+x^2*D2y-3*x*Dy=4*x^2','x')
ans =
C13 + C15*x^(3^(1/2) + 1) + C14*x^(1 - 3^(1/2)) - x^2 + (x^(1 - 3^(1/2))*x^(3^(1/2) + 1)*(3^(1/2)/3 - 1))/(3^(1/2) - 1) + (x^(1 - 3^(1/2))*x^(3^(1/2) + 1)*(3^(1/2)/3 + 1))/(3^(1/2) + 1)
2．求方程
3)0(',1)0(,'2")1(2===+y y xy y x
的解析解和数值解,并进行比较 解：解析解为
y=dsolve('(1+x^2)*D2y=2*x*Dy','y(0)=1','Dy(0)=3','x')
y =
x*(x^2 + 3) + 1 数值解：
设
12,'y y y y ==则原方程化为微分方程组
12222'2'1y y xy y x =⎧⎪⎨=⎪+⎩
定义函数m 文件fun7_2.m 如下 function f=fun7_2(x,y)
f=[y(2);2*x.*y(2)./(1+x.^2)]; 再用ode45求解
[x,y]=ode45('fun7_2',[0,5],[1,3]); plot(x,y(:,1),'ro') hold on
ezplot('x*(x^2 + 3) + 1',[0,5])
3．分别用ode45和ode15s 求解Van-del-Pol 方程
()⎪⎩
⎪⎨⎧===---1)0',0)0(0
)1(1000222x x x dt dx
x dt
x d 的数值解,并进行比较.
解：设1
2,'x x x x ==则原方程化为微分方程组
12
2
2121
''1000(1)x x x x x x =⎧⎨=-+⎩ 定义函数m 文件fun7_3.m 如下
function f=fun7_3(t,x)
f=[x(2);1000*(1-x(1).^2).*x(2)+x(1)];
再用ode45和ode15s 分别求解此方程，并绘图比较 clear;clf
[t1,x1]=ode45('fun7_3',[0,0.1],[0,1]); [t2,x2]=ode15s('fun7_3',[0,0.1],[0,1]); plot(t1,x1(:,1),'ro',t2,x2(:,1))
4．(单摆运动的近似解析解)当单摆初始角度0θ较小时,)(0θθ≤也较小,从θθ≈sin ,单摆运动微分方
程可近似写为
0)0(',)0(,"0===θθθθθmg ml
求此方程的解析解,并与练习3中的数值解进行比较.
解：用Matlab 命令求此方程的解析解，按练习3中的取值15)0(,8.9,1===θg l clear;close;
s=dsolve('D2y=9.8*y','y(0)=1*pi/180','Dy(0)=0','t')
s =
pi/(360*exp((7*5^(1/2)*t)/5)) + (pi*exp((7*5^(1/2)*t)/5))/360
练习3的中数值解 %M 文件fun7_4.m
function f=fun7_4(t,y)
f=[y(2), 9.8*sin(y(1))]'; %f 向量必须为一列向量 运行MATLAB 代码
[t,y]=ode45('fun7_4',[0,10],[1*pi/180,0]); s=eval(s);
plot(t,y(:,1),'ro');
xlabel('t'),ylabel('y1') hold on plot(t,s)
xlabel('t'),ylabel('y2')
clear;close;
s=dsolve('D2y=9.8*y','y(0)=1*pi/180','Dy(0)=0','t') [t,y]=ode45('fun7_4',[0,2],[1*pi/180,0]);
s=eval(s);
plot(t,y(:,1),'ro');
hold on
plot(t,s)
s =
由图象可知，区间改为[0,2]上时能观察出大概在区间[0,1.5]内二者能够较好吻合。

实验八平面图形的几何变换
【练习与思考】
1.将函数
2
x
y e-
=的图形向右平移3个单位且向上平移3个单位.
解：Matlab代码为
clear; close;
x=-2:0.1:2;y=exp(-x.^2);
x1=x+3; %图形向右平移3个单位;
y1=y+3;%图形向上平移3个单位; plot(x,y,x1,y,':',x1,y1,':','linewidth',3); axis([-3,6,-1,5])
xlabel('x'); ylabel('y');
grid on
2. 将函数2
x y e -=的图形在水平方向收缩一倍，在垂直方向放大一倍。

clear; close;
x=-2:0.1:2;y=exp(-x.^2);
x1=x/2; %图形在水平方向收缩一倍; y1=y*2;%图形在垂直方向放大一倍 plot(x,y,x1,y,'-.',x1,y1,':','linewidth',3); axis([-3,3,-1,3])
xlabel('x'); ylabel('y'); grid on
3. 将函数2
y x =的图形以原点为中心，顺时针旋转30度角. clear; close; x=-2:0.1:2; y=x.^2;
x1=x*cos(-pi/6)-y*sin(-pi/6); y1=x*sin(-pi/6)+y*cos(-pi/6); plot(x,y,x1,y1,'r:','linewidth',3);
legend('原图','顺时针旋转30度角后的图') xlabel('x'); ylabel('y'); grid on
4. 已知函数
2,2,02y x x x =-≤≤，试扩展函数的定义域，使之成为以2周期的偶函数，并画出函数
在[-8，8]上的图形。

若要把函数延拓成以4为周期的奇函数呢?
解：延拓成以2周期的偶函数，画出函数在[-8，8]上的图形的Matlab 代码为： clear; close;
x=0:0.1:2;y=2*x-x.^2; x1=-x;
y1=-2*x1-x1.^2; x2=x+2;y2=y; x3=-x-2;y3=y1; x4=x2+2;y4=y2;
x5=x3-2;y5=y3; x6=x4+2;y6=y4; x7=x5-2;y7=y5;
plot(x,y,x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7); xlabel('x'); ylabel('y'); 把函数延拓成以4为周期的奇函数，画出函数在[-8，8]上的图形的Matlab 代码为： clear; close;
x=0:0.1:2;y=2*x-x.^2; x1=-x;
y1=2*x1+x1.^2; x2=x-4;y2=y; x3=x1+4;y3=y1; x4=x1-4;y4=y1; x5=x2+8;y5=y2;
x6=x2-4;y6=y2; x7=x3+4;y7=y3;
plot(x,y,x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7);
%legend('x-y','x1-y1','x2-y2','x3-y3','x4-y4','x5-y5','x6-y6','x7-y7') xlabel('x'); ylabel('y'); grid on
5. 做怎样的变换才能使函数图形绕给定的点(,)a b 转动?这个变换可以分解成3个基本变换:平移量为
(,)a b --的平移变换1T ,旋转角度为α的旋转变换2T ,1T 的逆变换11T -.求出变换矩阵,写出与变换相应的方程,
并对具体的函数图形进行变换. a)
sin ,(0,2)y x x π=∈ (2) sin ,cos ,(0,2)x a t y b t t π==∈
解：a)
clear; close; a=pi;b=0;
alpha=60*pi/180;
T1=[1 0 -a;0 1 -b;0 0 1];
T2=[cos(alpha) -sin(alpha) 0; sin(alpha) cos(alpha) 0; 0 0 1];
T=inv(T1)*T2*T1; %inv 求矩阵的逆 x=0:0.1*pi:2*pi;y=sin(x);
x1=T(1,1)*x+T(1,2)*y+T(1,3); y1=T(2,1)*x+T(2,2)*y+T(2,3);
plot(x,y,x1,y1,a,b,'.','markersize',35);
xlabel('x'); ylabel('y');
text(a,b,'\leftarrow(a,b)','fontsize',30) grid on。