数学实验课后习题解答

实验一 曲线绘图 【练习与思考】

画出下列常见曲线的图形。 以直角坐标方程表示的曲线：

1. clear;

x=-2:0.1:2; y=x.^3; plot(x,y)

3

y x立方曲线

y=3*a*t.^3./(1+t.^2);

plot(x,y)

8. 摆线x a(t sint),yclear;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. 内摆线（星形线）

b(1 cost)

2. 立方抛物线clear;

y=-2:0.1:2; x=y.^3; plot(x,y) grid on

y x

x acost,y asint(x y a)

2

33

232323

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

x

y e高斯曲线

clear; a=1;

t=0:pi/50:2*pi; x=a*cos(t).^3; y=a*sin(t).^3; plot(x,y)

10. 圆的渐伸线（渐开线）

x a(cost tsint),y a(sint tcost)

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. 空间螺线

x t,y t(y x)

32

23

x acost,y bsint,z ct

t2,y t3(y2 x3)

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. 阿基米德线r a ,rclear; a=1;

phy=0:pi/50:6*pi; rho=a*phy;

polar(phy,rho,'r-*') 13. 对数螺线r eclear; a=0.1;

phy=0:pi/50:6*pi; rho=exp(a*phy); polar(phy,rho) 14. 双纽线

a

5. 半立方抛物线xclear;

t=-3:0.05:3; x=t.^2;y=t.^3; plot(x,y) %axis equal grid on

6. 迪卡尔曲线

x

3at3at33

,y (x y 3axy 0)22

1 t1 t

2

0

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. 蔓叶线

at2at3x32

x ,y (y )

1 t21 t2a x

clear;

a=3;t=-6:0.1:6;

x=3*a*t.^2./(1+t.^2);

r2 a2cos2 ((x2 y2)2 a2(x2 y2))

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. 双纽线

16. 四叶玫瑰线r clear;close a=1;

phy=0:pi/50:2*pi; rho=a*sin(2*phy); polar(phy,rho)

17. 三叶玫瑰线r clear;close a=1;

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

asin2 ,r 0

asin3 ,r 0

r2 a2sin2 ((x2 y2)2 2a2xy)

clear; a=1;

phy=0:pi/50:pi/2;

rho=a*sqrt(sin(2*phy)); polar(phy,rho) hold on

polar(phy,-rho)

a=1;

phy=0:pi/50:2*pi; rho=a*cos(3*phy); polar(phy,rho)

acos3 ,r 0

实验二 极限与导数

【练习与思考】

1． 求下列各极限

（1）lim(1

n

1n

) （2）limn3 3n

n n

（3）lim(

n

n 2 2n 1 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）lim(

x 1

21

) （5）limxcot2x （6）

x 0x2 1x 1

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． 考虑函数

lim(x2 3x 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）lim(cos

f(x) 3x2sin(x3), 2 x 2

作出图形，并说出大致单调区间；使用diff求f'(x)，并求f(x)确切的单调区间。

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)

mx11

) （8）lim( x) （9）

x x 1xxe 1 x 1lim x 0x

clear;

用fzero函数找f'(x)的零点，即原函数f(x)的驻点 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) 求

f'(x)所有零点（即f(x)的驻点）； (iii) 求出驻点处f(x)的二阶导数值;

(iv) 用fmin求各极值点的确切位置; (v) 局部极值点与f'(x),f"(x)有何关系?

(1) f(x) x2

sin(x2

x 2),x [ 2,2] (2) f(x) 3x5 20x3 10,x [ 3,3] (3)

f(x) x3

x2

x 2,x [0,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函数找f'(x)的零点，即原函数f(x)的驻点 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命令观测函数

y f(x)的Maclaurin展开式的前

几项, 然后在同一坐标系里作出函数y f(x)和它的Taylor展开

式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向

y f(x)的图形的逼近的情况 (1) f(x) arcsinx

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) f(x) arctanx 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) f(x) ex2

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)

f(x) sin2x

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 x

(5)

f(x)

e1 x

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 + (163*x^5)/60 + (65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1

(6) f(x) ln(x x2

) 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. 求公式 1 2k

,k 1,n 1n

2k m2, ,中的数kmk,k 4,5,6,7,8的值.

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. 利用公式 1

e来计算e的近似值。精确到小数点后n 0n!

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)

1n 1n2 n

3

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)

1n 1n2

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) sin1

n 1

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)

lnn

3 n 1

nclear;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)

n!

n 1

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)

1nn 3(lnn)

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

n 1

nlnn 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

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)

plot(s,'.')

(8) ( 1)nnn 1n2 1

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验证

xsinx2

dx， dx

1 cosx，

dx

ex 1，

arcsinxdx， sec3

xdx

解：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)

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

sinxx

2

x

0xdx,

1

xdx, esin(2x)dx, 1

e x2

dx

解：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

x23．(椭圆的周长) 用定积分的方法计算椭圆y2

9 4

1的周长

解：椭圆的参数方程为

x 3cost

y 2sint由参数曲线的弧长公式得

s

2

2

2 0

Matlab代码为

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

s =

15.865

【练习与思考】

4．（二重积分）计算数值积分

(1 x y)dxdy

x2 y2 2y

解：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

1/3

1

x

cosxdx，会出现什么问题？分析原因，并求出正确的

解。

解：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．（假收敛现象）考虑积分I(k) k

sinxdx，

（1）用解析法求I(k);

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求I(4),I(6)和I(8),发现什么问题？

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法数值积分程序.计算公式为

I Sh

n

3

(f1 4f2 2f3 4f4 2fn 1 4fn

其中n为偶数,h

b a

n

,fi f(a (i 1)h),i 1,2, ,n 1. 解：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．(广义积分)计算广义积分

exp( x2)

1tan(x)1, sinx 1 x4dx, 0xdx0 x2dx并验证公式

x2

exp( )sinx 2

dx 1,

xdx

2

.

解：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 x,y 30范围内的图形，并画出相应的等高

线。 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) 椭球面

x 3cosusinv,y 2cosucosv,z sinu

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)

椭圆抛物面

x 3usinv,y 2ucosv,z 4u2

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) 单叶双曲面

x 3secusinv,y 2secucosv,z 4tanu

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) 2d)

u,y v,z

u v2

双曲抛物面x3

clear

u=-3:0.1:3;

[U,V]=meshgrid(u); x=U; y=V;

z=(U.^2-V.^2)/3; mesh(x,y,z)

clc;clear; x lnusinv,y lnucosv,z u

x1=[100 100 100 100 200 200 200 200 300 300 300 300 400 400 400

clear; 400]; u=1:0.1:5; x2=[100 200 300 400 100 200 300 400 100 200 300 400 100 200 300 v=0:pi/30:2*pi; 400]; [U,V]=meshgrid(u,v); y=[636 698 680 662 697 712 674 626 624 630 598 552 478 478 412 x=log(U).*sin(V); 334]'; y=log(U).*cos(V); x=[x1',x2']; z=U; x0=[1 1 1 1 1]; mesh(x,y,z) beta=lsqcurvefit('heigh',x0,x,y) axis equal %绘图：

f) 圆锥面x usinv,y ucosv,z u a1=100:5:400;

a2=a1;

[xx1,xx2]=meshgrid(a1,a2); clear;

Z=beta(1)+beta(2)*xx1+beta(3)*xx2+beta(4)*xx1.^2+beta(5)*xx2.^u=-5:0.1:5;

2; v=0:pi/30:2*pi;

mesh(xx1,xx2,Z) [U,V]=meshgrid(u,v);

x=(U).*sin(V);

y=(U).*cos(V);

Local minimum possible. z=U;

mesh(x,y,z)

lsqcurvefit stopped because the final change axis equal

in the sum of squares relative to

its initial value is less than the default g) 环面

value of the function tolerance.

x (3 0.4cosu)cosv,y (3 0.4cosu)sinv,z 0.4sinv clear; u=0:pi/30:2*pi; beta = v=u; Columns 1 through 5

538.4375 1.4901 0.6189 -0.0046 [U,V]=meshgrid(u,v);

-0.0017 x=(3+0.4*cos(U)).*cos(V);

y=(3+0.4*cos(U)).*sin(V); z=0.4*sin(V); contour(xx1,xx2,Z,30),colorbar mesh(x,y,z) %计算最高点及高程 x0=[100,100];

options=optimset('largescale','off'); h) 正螺面x usinv,y ucosv,z 4v

%设置下界

lb=[0,0]; clear;

%无上界 u=0:pi/30:pi;

ub=[]; v=0:pi/30:10*pi;

[x,fval]=fmincon('height',x0,[],[],[],[],lb,ub,[],options) [U,V]=meshgrid(u,v);

x=U.*sin(V);

Warning: Options LargeScale = 'off' and y=U.*cos(V);

Algorithm = 'trust-region-reflective' z=4*V;

conflict. mesh(x,y,z)

Ignoring Algorithm and running active-set colorbar

algorithm. To run trust-region-reflective,

set 3. 在一丘陵地带测量搞程,x和y方向每隔100米测一个点,得

LargeScale = 'on'. To run active-set without

高程见表5-2,试拟合一曲面,确定合适的模型,并由此找出最高点和

this warning, use Algorithm = 'active-set'.

该点的高程.

> In fmincon at 445 e)

旋转面

were satisfied to within the default value of the constraint tolerance. No active inequalities. x =

161.9676 182.0320 fval =

-715.4403

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.*xx1+beta(6)*xx2.^2; end

%height.m

function y=height(x)

y=-(434.0000+1.9079*x(1)+1.0366*x(2)-0.0017*x(1).^2-0.0046*x(2).*x(1)-0.0017*x(2).^2); end

实验六 多元函数的极值

【练习与思考】 1.求z

x4 y4 4xy 1的极值，并对图形进行观测。

解：Maltab代码为 syms x y;

z=x^4+y^4-4*x*y+1; dzx=diff(z,x); dzy=diff(z,y);

[x,y]=solve(dzx,dzy,x,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.求函数

f x,y x2 2y2在圆周x2 y2 1的最大值和最

小值。

解：构造Lagrange函数

L(x,y) x2 2y2 (x2 y2 1)

求Lagrange函数的自由极值.先求L关于x,y, 的一阶偏导数，

再解正规方程可得所求的极值点，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）为函数的条件极

值点，经判断函数f x,y x2 2y2在（1,0）、（-

1,0）取得极小值，在（0,1）、（0，-1）取得极大值。

3.在球面x

2

y2 z2 1求出与点(3,1,-1)距离最近和最远点。

解：设球面上的点为(x,y,z),则此点与点(3,1,-1)的距离为

d x,y,z (x 3)2 (y 1)2 (z 1)2

且(x,y,z)

满足x

2

y2 z2 1；构造Lagrange函数

L(x,y,z, ) (x 3)2 (y 1)2 (z 1)2 (x2 y2

求Lagrange函数的自由极值.先求L关于x,

y,z, 的一阶偏导

数，再解正规方程可得所求的极值点，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]

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]

经判断可知，函数

f(x,y,z) x 2y 3z

在平面

x y z 1与柱面x

2 y2 1的交线上的最大值为3 。

22

5.求函数z x y在三条直线x 1,y 1,x y 1所围

区域上的最大值和最小值。

解：显然此函数的驻点为(0,0)不在此区域内，因此该函数的最大值和最小值点应在三条边界上，下面分别求此函数在这三条边界上的最大值和最小值，Matlab代码如下

（1）求函数在直线边界x=1,0 y 1上的最大值和最小值

将x=1代入原函数，则二元函数变为一元函数

z=1+y2,0 y 1

最大值点为y=1，最大值为2，最小值点为y=0，最小值为1； （2）求函数在直线边界y=1,0 x 1上的最大值和最小值

将x=1代入原函数，则二元函数变为一元函数

z=1+x2,0 x 1

最大值点为x=1，最大值为2，最小值点为x=0，最小值为1； (3)求函数在直线边界x+y=1,0 x 1上的最大值和最小值

将y=1-x代入原函数，则二元函数变为一元函数

z=(1-x)2+x2,0 x 1

用Matlab命令求此函数的最大和最小值点 先求驻点 clear; syms x y2 1)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。

、得到条件极值点为( ，经判断，球面

111111

x2 y2 z2 1上与点(3,1,-1)

距离最近的点为,

，最远的点111111(。

(

4.求函数面x

2

f(x,y,z) x 2y 3z在平面x y z 1与柱

y2 1的交线上的最大值。

解：构造Lagrange函数

L(x,y,z, ) x 2y 3z 1(x y z 1) 2(x2

求Lagrange函数的自由极值.先求L关于x,y,z, 1, 2的一阶偏导数，再解正规方程可得所求的极值点，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'

实验七 常微分方程

【练习与思考】

1．求下列微分方程的解析解

a) 一阶线性方程y' x3

y 2 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) 贝努利方程y' xy2 y 0 dsolve('Dy-x*y^2-y=0','x')

ans =

0 exp(x)/(C4 - exp(x)*(x - 1))

c) 高阶线性齐次方程y"' y" 3y' 2y 0 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) 高阶线性非齐次方程y" 3y' 2y 3sinx 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) 欧拉方程x

3

y"' x2y" 3xy' 4x2

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．求方程

(1 x2)y" 2xy',y(0) 1,y'(0) 3

的解析解和数值解,并进行比较 解：解析解为

y=dsolve('(1+x^2)*D2y=2*x*Dy','y(0)=1','Dy(0)=3','x')

y =

x*(x^2 + 3) + 1 数值解：

设

y1 y,y2 y'则原方程化为微分方程组

y1' y2 y2xy2

2' 1 x2定义函数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方程

d2x2dx

2 1000(1 x) dt

dt x 0

x(0) 0,x' 0) 1 的数值解,并进行比较.

解：设x1

x,x2 x'则原方程化为微分方程组

x1' x2

' 1000(1 x x x2

21)x21

定义函数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

,单摆运动微分方程可近似写为ml " mg , (0) 0, '(0) 0

求此方程的解析解,并与练习3中的数值解进行比较.

解：用Matlab命令求此方程的解析解，按练习3中的取值

l 1,g 9.8, (0) 15

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 =

pi/(360*exp((7*5^(1/2)*t)/5)) + (pi*exp((7*5^(1/2)*t)/5))/360

由图象可知，区间改为[0,2]上时能观察出大概在区间[0,1.5]内二者能够较好吻合。

实验八 平面图形的几何变换

【练习与思考】

1. 将函数

y e x

2

的图形向右平移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. 将函数

y e x

2

的图形在水平方向收缩一倍，在垂直方向

放大一倍。

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. 将函数y x2

的图形以原点为中心，顺时针旋转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. 已知函数

y 2x x2,,0 x 2，试扩展函数的定义

域，使之成为以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)的平移变换T1,旋转角度为 的旋转变换T 1

2,T1的逆变换T1.求出变换矩阵,写出与变换相应的方程,并对具体的函数图形进行变换. a)

y sinx,x (0,2 ) (2)

x asint,y bcost,t (0,2 )

解：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 b)

clear; close; a=1;b=2;

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]; y3=sin(2*x+cos(x+sin(x))); plot(x,y3) grid on

2. 求出练习1中5个函数的导函数，画图研究这些函的周期性，并从理论上证明周期函数的导函数仍为周期函数。 （1）sin2x close;clc;clear syms x;

y=sin(x+sin(2*x)); T=inv(T1)*T2*T1; %inv求矩阵的逆 t=0:0.1*pi:2*pi; x=a*sin(t); y=b*cos(t);

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(a,b,'o',x,y,x1,y1);

xlabel('x'); ylabel('y'); grid on

实验十 周期函数

【练习与思考】

1. 画图研究下列函的周期性，并从理论上证明：（

1

）

sin2x

（cosx 2cos2x 3cos3x 4cos4x （3）sin(x cos(x sin(x))) （4）sin(x sin(x sin(x))) （5）sin(2x cos(x sin(x)) 解：

clear;close;

x=-3*pi:0.05*pi:3*pi; y1=sin(2*x); plot(x,y1) grid on

y2=cos(x)+2*cos(2*x)+3*cos(3*x)+4*cos(4*x); plot(x,y2)

axis([-12,12,-7,11]) grid on

（3）sin(x cos(x sin(x))) clear;close;

x=-4*pi:pi/50:4*pi;

y3=sin(x+cos(x+sin(x))); plot(x,y3) grid on

（4）sin(x sin(x sin(x))) clear;close;

x=-4*pi:pi/50:4*pi;

y3=sin(x+sin(x+sin(x))); plot(x,y3) grid on

（5）sin(2x cos(x sin(x)) clear;close;

x=-4*pi:pi/50:4*pi;

y1=diff(y);

x=-2*pi:pi/30:2*pi; y=subs(y1,x); plot(x,y)

grid on

（2）cosx 2cos2x 3cos3x 4cos4x close;clc;clear syms x;

y=cos(x)+2*cos(2*x)+3*cos(3*x)+4*cos(4*x); y1=diff(y);

x=-2*pi:pi/30:2*pi; 2

）

y=subs(y1,x); plot(x,y) grid on

（3）sin(x cos(x sin(x))) close;clc;clear syms x;

y=sin(x+cos(x+sin(x))); y1=diff(y);

x=-2*pi:pi/30:2*pi; y=subs(y1,x); plot(x,y)

grid on

（4）sin(x sin(x sin(x))) close;clc;clear syms x;

y=sin(x+sin(x+sin(x))); y1=diff(y);

x=-4*pi:pi/30:4*pi; y=subs(y1,x); plot(x,y)

grid on

（5）sin(2x cos(x sin(x)) close;clc;clear syms x;

y=sin(2*x+cos(x+sin(x))); y1=diff(y);

x=-4*pi:pi/30:4*pi; y=subs(y1,x); plot(x,y)

grid on

本文来源：https://www.bwwdw.com/article/ga1m.html