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

2024年1月17日发(作者：西藏小学内地班数学试卷)

数学实验课后习题解答

配套教材：王向东 戎海武 文翰 编著数学实验

王汝军编写

实验一 曲线绘图

【练习与思考】

画出下列常见曲线的图形。

以直角坐标方程表示的曲线：

1.

clear;

x=-2:0.1:2;

y=x.^3;

plot(x,y)

86420-2-4-6-8-23yx立方曲线

-1.5-1-0.500.511.52

2. 立方抛物线clear;

y=-2:0.1:2;

x=y.^3;

plot(x,y)

grid on

21.510.50-0.5-1-1.5-2-8

y3x

-6-4-202468

3. 高斯曲线clear;

x=-3:0.1:3;

y=exp(-x.^2);

plot(x,y);

grid on

%axis equal

yex2

10.90.80.70.60.50.40.30.20.10-3-2-10123

以参数方程表示的曲线

4. 奈尔抛物线clear;

t=-3:0.05:3;

x=t.^3;y=t.^2;

plot(x,y)

axis equal

grid on

2520151050-5-10-15-25-20-15-10-5

23xt3,yt2(yx)

5. 半立方抛物线xclear;

t=-3:0.05:3;

x=t.^2;y=t.^3;

plot(x,y)

%axis equal

grid on

3020

t2,yt3(y2x3)

100-10-20-3

6.

3at3at233,y(xy3axy0) 迪卡尔曲线x221t1tclear;

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)

9876543210-5-4-3-2-1012345

7.

at2at3x32蔓叶线x,y(y)

221t1taxclear;

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)

6040200-20-40-6

8. 摆线xa(tclear;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

8642

sint),yb(1cost)

0-2-4-6618

9. 内摆线（星形线）xclear;

acos3t,yasin3t(xya)

232323

a=1;

t=0:pi/50:2*pi;

x=a*cos(t).^3;

y=a*sin(t).^3;

plot(x,y)

10.80.60.40.20-0.2-0.4-0.6-0.8-1-1-0.8-0.6-0.4-0.200.20.40.60.81

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

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

20151050-5-10-15-20-20

a(costtsint),ya(sinttcost)

-15-10-5051015

11. 空间螺线xclear

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

acost,ybsint,zct

2-1-2-4-2204

以极坐标方程表示的曲线：

12. 阿基米德线r

a,r0

clear;

a=1;

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

rho=a*phy;

polar(phy,rho,\'r-*\')

90 20120 15150 10 5

13. 对数螺线reclear;

a=0.1;

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

rho=exp(a*phy);

polar(phy,rho)

120 6150 4 218030

a

90 8670300

14. 双纽线racos2clear;

a=1;

phy=-pi/4:pi/50:pi/4;

rho=a*sqrt(cos(2*phy));

polar(phy,rho)

22

((x2y2)2a2(x2y2))

hold on

polar(phy,-rho)

90120 0.8 0.6150 0.4 0.2180030 160300

15. 双纽线rasin2clear;

a=1;

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

rho=a*sqrt(sin(2*phy));

polar(phy,rho)

hold on

polar(phy,-rho)

90120 0.8 0.6150 0.4 0.2180030 160

22((x2y2)22a2xy)

210

16. 四叶玫瑰线rclear;close

a=1;

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

rho=a*sin(2*phy);

polar(phy,rho)

90120

asin2,r0

160 0.8 0.6150 0.4 0.2340270300

17. 三叶玫瑰线rclear;close

a=1;

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

rho=a*sin(3*phy);

polar(phy,rho)

asin3,r0

90120 160 0.8 0.6150 0.4 0.2340270300

18. 三叶玫瑰线rclear;close

a=1;

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

rho=a*cos(3*phy);

polar(phy,rho)

acos3,r0

90120 160 0.8 0.6150 0.4 0.2340270300

【练习与思考】

1． 求下列各极限

（1）lim(1n

实验二 极限与导数

1n) （2）limnn33nn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1212limxcot2x （5） （6）)lim(x3xx)

2x0xx1x1

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

31x1mx11（7）lim(cos) （8）lim(x

) （9）limx0xx1xxxe1clear;

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

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];

3 x2 sin(x3)1050-5-10-2-1.5-1-0.50x0.511.52

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)

6 x sin(x3) + 9 x4 cos(x3)120-20-40-60-2-1.5-1-0.5

用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) 作出图形，观测所有的局部极大、局部极小和全局最大、全局最小值点的粗略位置；

0x0.511.52f\'(x)所有零点（即f(x)的驻点）；

(iii) 求出驻点处f(x)的二阶导数值;

(iI) 求(iv) 用fmin求各极值点的确切位置;

(v) 局部极值点与(1)

(2)

(3)

f\'(x),f\"(x)有何关系?

f(x)x2sin(x2x2),x[2,2]

f(x)3x520x310,x[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)

x2 sin(x2 - x - 2)210-1-2-3-2-1.5-1-0.50x0.511.52

局部极大值点为：-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)

2 x sin(x2 - x - 2) + x2 cos(x2 - x - 2) (2 x - 1)20151050-5-2-1.5-1-0.5

用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

0x0.511.52

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

21.510.50-0.5-1-1.5-2-1-0.8-0.6-0.4-0.200.20.40.60.81

(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

10.80.60.40.20-0.2-0.4-0.6-0.8-1-1-0.8-0.6-0.4-0.200.20.40.60.81

x2

(3)

f(x)e

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

2.82.62.42.221.81.61.41.21-1-0.8-0.6-0.4-0.200.20.40.60.81

2

(4)

f(x)sinx

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

50-5-10-15-20-25-4-3-2-101234

(5)

exf(x)1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)/ + (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

21.510.50-0.5-1-1.5-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.10

2

(6)

f(x)ln(x1x)

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

10.80.60.40.20-0.2-0.4-0.6-0.8-1-1-0.8-0.6-0.4-0.200.20.40.60.81



12k,k1,2,,中的数mk,k4,5,6,7,8的值. 2. 求公式2kmnn1kk=[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的近似值。精确到小数点后100位,这时应计算到这个无穷级数的前多少n!n0项?请说明你的理由.

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

n =

70

精确到小数点后100位,这时应计算到这个无穷级数的前71项,理由是误差小于10的负100次方，需要最后一项小于10的负100次方，由上述循环知n=70时最后一项小于10的负100次方，故应计算到这个无穷级数的前71项.

4. 用练习3中所用观测法判断下列级数的敛散性

(1)

123n1n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

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,\'.\')

0.660.640.620.60.580.560.540.520.5

(2)

1nn1n2

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

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,\'.\')

0.70.680.660.640.620.60.580.560.540.520.5

(3)

sinn11

nclear;clc;

epsl=0.01;

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)

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,\'.\')

200-20-40-60-80-100-126

发散

(4)

lnn

3n1n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(\'收敛\')

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,\'.\')

0.20.180.160.140.120.10.080.060.040.5

(5)

n!

nn1nclear;close

syms n

s=[];he=0;

for k=1:100

he=he+factorial(k)/k^k;

s(k)=he;

end

plot(s,\'.\')

1.91.81.71.61.51.41.31.21.11

(6)

1nn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

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,\'.\')

1.41.210.80.60.40.26

(7)

1

nlnnn1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)

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,\'.\')

32.521.510.50250300

(1)nn (8)

2

n1n1clear;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)

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,\'.\')

0.350.30.250.20.150.10.200250300

实验四 积分

【练习与思考】

1．(不定积分)用int计算下列不定积分，并用diff验证

2xsinxdx，dx1cosx，dx3arcsinxdxsec，，xdx

ex1解：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计算下列定积分

sinx0xdx,110xdx,esin(2x)dx,exdx

00x2x12解：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

x2y23．(椭圆的周长) 用定积分的方法计算椭圆1的周长

94x3cost解：椭圆的参数方程为

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

s20x\'(t)2y\'(t)2dt209sin2t4cos2tdt205sin2t4dt

Matlab代码为

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

s =

15.865

4．（二重积分）计算数值积分2xy22y(1xy)dxdy

解：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/3xcosxdx，会出现什么问题？分析原因，并求出正确的11

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

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

h(f14f22f34f42fn14fnfn1)

3ba其中n为偶数,h,fif(a(i1)h),i1,2,,n1.

nISn解：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．(广义积分)计算广义积分

1sinx1tan(x)exp(x2)1x4dx,0xdx,01x2dx

并验证公式

x2exp()2dx1,20sinxdxx2.

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

实验五 二元函数的图形

【练习与思考】

1. 画出空间曲线z10sinx2y21xy22在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)

clear;

椭球面x3cosusinv,y2cosucosv,zsinu

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) 椭圆抛物面xclear;

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

3usinv,y2ucosv,z4u2

c) 单叶双曲面x3secusinv,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)

y2secucosv,z4tanu

d)

u2v2双曲抛物面xu,yv,z

3clear

u=-3:0.1:3;

[U,V]=meshgrid(u);

x=U;

y=V;

z=(U.^2-V.^2)/3;

mesh(x,y,z)

e) 旋转面xlnusinv,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

ylnucosv,zu

f)

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

圆锥面xusinv,

yucosv,zu

g) 环面x(30.4cosu)cosv,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)

y(30.4cosu)sinv,z0.4sinv

h)

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

正螺面xusinv,yucosv,z4v

3. 在一丘陵地带测量搞程,x和y方向每隔100米测一个点,得高程见表5-2,试拟合一曲面,确定合适的模型,并由此找出最高点和该点的高程.

表5-2 高程数据

y x

100 200 300 400

100

200

300

400

636

698

680

662

697

712

674

626

624

630

598

552

478

478

412

334

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 =

-715.4403

400

700

100350400

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.求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);

s=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)

21.81.61.41.2110.50-0.5-1-1-0.50.501

可得点（1,0）、（0,1）（-1,0）、（0，-1）为函数的条件极值点，经判断函数（1,0）、（-1,0）取得极小值，在（0,1）、（0，-1）取得极大值。

3.在球面x解数

：2fx,yx22y2在y2z21求出与点(3,1,-1)距离最近和最远点。

球面上的点为(x,y,z),则此点与点(3,1,-1)的距离为Lagrange函设dx,y,z(x3)2(y1)2(z1)2且(x,y,z)满足x2y2z21；构造L(x,y,z,)(x3)2(y1)2(z1)2(x2y2z21)

求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]

得到条件极值点为(311,,)、(,,)，经判断，球面13111111,,)，最远的点x2y2z21上与点(3,1,-1)距离最近的点为(11(,,)。

111111

4.求函数f(x,y,z)x2y3z在平面xyz1与柱面x2y21的交线上的最大值。

解：构造Lagrange函数

L(x,y,z,)x2y3z1(xyz1)2(x2y21)

求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\'

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]

经判断可知，函数大值为3f(x,y,z)x2y3z在平面xyz1与柱面x2y21的交线上的最29。

225.求函数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

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) 一阶线性方程y\'xy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) 贝努利方程3y\'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

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

3y\"\'x2y\"3xy\'4x2

(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2xy2y\'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])

x (x2 + 3) + 0.511.522.5x33.544.55

3．分别用ode45和ode15s求解Van-del-Pol方程

d2x2dx21000(1x)x0

dtdtx(0)0,x\'0)1的数值解,并进行比较.

解：设x1x,x2x\'则原方程化为微分方程组

x1\'x2

2x\'1000(1x)xx2121定义函数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))

1.81.61.41.210.80.60.40.2000.010.020.030.040.050.060.070.080.090.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\')

3.5x 101132.52y21.510.50012345t678910显然，在这个区间内，二者差别较大，下面改为较小区间

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

54.543.532.521.510.5000.20.40.60.811.21.41.61.82

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

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

【练习与思考】

1. 将函数yex2的图形向右平移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

5432y10-1-3-2-101x23456

x2

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

32.521.510.50-0.5-1-3y-2-10x123

2

3. 将函数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.543.532.521.510.50-0.5

-2-101x23y

原图顺时针旋转30度角后的图4

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\');

10.90.80.70.60.50.40.30.20.10-8-6-4-20x2468y

把函数延拓成以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

10.80.60.40.20-0.2-0.4-0.6-0.8-1-8-6-4-20x2468y

5. 做怎样的变换才能使函数图形绕给定的点(a,b)转动?这个变换可以分解成3个基本变换:平移量为(a,b)的平移变换T1,旋转角度为的旋转变换T2,T1的逆变换T11.求出变换矩阵,写出与变换相应的方程,并对具体的函数图形进行变换.

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

3210(a,b)y-1-2-30123x4567

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];

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

21.510.50-0.5-1-1.5-2-1y012x345

的近似计算 实验九

【练习与思考】

1. 利用勾股定理推导在割圆术中给出的公式

2x62n124x6,S62n!32nx62n，

S2n2S2nSn。

2n解：略

2. 利用韦达公式，构造出一种迭代算法来计算的近似值，并进行实际计算，评价算法效果。

解：韦达(VieTa)公式

22设a12222222222222

2,an2an1,(n2,3,4,),s02,sn2sn1,(n1,2,3,),则limsnnan于是得到一种迭代算法，实际计算的Matlab代码为

clear;

a=sqrt(2);

s=2;

n=15;

for k=1:n

s=2*s/a;

a=sqrt(2+a);

pai=vpa(sym(s),20)

error=vpa(sym(pi)-sym(s),20)

end

pai =

2.82842776

error =

0.360314086

pai =

3.7178101

error =

0.075428318

pai =

3.80519693

error =

0.741269122

pai =

3.13654849

error =

0.854366431

pai =

3.1425117

error =

0.

pai =

3.14127725

error =

0.

pai =

3.1415138

error =

0.492339280772

pai =

3.14157294

error =

0.7

pai =

3.71595707

error =

0.633667782247

pai =

3.11997443

error =

0.5934941751198

pai =

3.01175726

error =

0.67566585639321

pai =

3.48724538

error =

0.927

pai =

3.85631725

error =

0.23

pai =

3.69856301

error =

0.866

pai =

3.65911335

error =

0.2649

3. 利用欧拉公式

2812n0(2n1)

计算出的前4位小数。

解:Matlab代码为

tic

clear;clc

syms n

k=1;

flag=1;

s=1;

while flag>0.00001

s=s+1/(2*k+1)^2;

pai=sqrt(s*8);

k=k+1;

flag=pi-pai;

end

k-1

vpa(pai,5)

toc

ans =

31831

ans =

3.1416

Elapsed time is 0.097453 seconds.

4．使用高斯公式

48arctan和斯托梅公式

11132arctan20arctan4arctan57239

24arctan8arctan计算出的前60位小数。

5．用改进后的公式

1的前n1264032032(1)n(6n)!13591409545140134n，

33n(n!)(3n)!640320n01,2,3,,15项计算的值，并列出计算误差。

6．使用Bailey的迭代算法

y05(52),cn(25yn1)2,dn5yn11,3endn((7cn)23dn7cn),yn25yn1(1dn525en2)en2,

21n1yn152a0,anyn1an15[yn1(yn12yn15)].22

迭代计算4次，计算的近似值，观测近似效果。并用迭代误差估计公式代次数和迭代误差界的数据表（n20）。

实验十 周期函数

【练习与思考】

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

（1）sin2x （2）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

10.80.60.40.20-0.2-0.4-0.6-0.8-1-10-8-6-4-20246810n1165ne5计算出迭an

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

1086420-2-4-6-10-50510（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

10.80.60.40.20-0.2-0.4-0.6-0.8-1-15-10-5051015

（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

10.80.60.40.20-0.2-0.4-0.6-0.8-1-15-10-5051015

（5）sin(2xcos(xsin(x))

clear;close;

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

y3=sin(2*x+cos(x+sin(x)));

plot(x,y3)

grid on

10.80.60.40.20-0.2-0.4-0.6-0.8-1-15-10-5051015

2.

求出练习1中5个函数的导函数，画图研究这些函的周期性，并从理论上证明周期函数的导函数仍