数学实验课后习题解答

数学实验课后习题解答王汝军编写实验一曲线绘图【练习与思考】画出下列常见曲线的图形。以直角坐标方程表示的曲线：.立方曲线y=X3clear;x=-2:0.1:2;y=x∙^3;plot(x,y)2.立方抛物线，=3χclear;y=-2:0.1:2;χ=y.^3;plot(x,y)gridon.高斯曲线y=eT2clear;x=-3:0.1:3;y=exp(-x,^2);plot(x,y);

gridon%axisequal以参数方程表示的曲线2.奈尔抛物线x=t3,y=t2(y=x3)clear;t=-3:0.05:3;x=t.^3;y=t.^2;plot(x,y)axisequalgridon.半立方抛物线x=t2,y=t3(y2=x3)clear;t=-3:0.05:3;x=t.^2;y=t.^3;plot(x,y)%axisequalgridon.迪卡尔曲线X=H,y=把(X3+y3-3Oxy=0)1+12 1+12clear;a=3;t=-6:0.1:6;x=3*a*t./(1+t.^2);y=3*a*t∙^2∙Z(1+t^2);plot(x,y).蔓叶线X=旦,y=士(y2=旦)1+2 1+2a-Xclear;a=3;t=-6:0.1:6;x=3*a*t.^2.Z(1+t.^2);y=3*a*t.^3.Z(1+t.^2);plot(x,y).摆线X=a(t-sint),y=b(1-cost)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);axisequalgridon222X=222X=acos31,y=asin31(X3+y3=a3)clear;a=1;t=0:pi/50:2*pi;x=a*cos(t).^3;y=a*sin(t).^3;plot(x,y)tcost).圆的渐伸线(渐开线)X=a(CoSt+1Sint),y=atcost)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)gridon.空间螺线x=acost,y=bsint,Z=Ctcleara=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)gridon

以极坐标方程表示的曲线：.阿基米德线r=aφ,r≥0clear;a=1;phy=0:pi/50:6*pi;rho=a*phy;polar(phy,rho,'r-*').对数螺线r=eaφclear;a=0.1;phy=0:pi/50:6*pi;rho=exp(a*phy);polar(phy,rho).双纽线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)holdonpolar(phy,-rho).双纽线r2=a2sin2φ((X2+y2)2=2a2χy)clear;a=1;phy=0:pi/50:pi/2;rho=a*sqrt(sin(2*phy));polar(phy,rho)holdonpolar(phy,-rho).四叶玫瑰线r=asin2φ,r≥0clear;closea=1;phy=0:pi/50:2*pi;rho=a*sin(2*phy);polar(phy,rho)

.三叶玫瑰线r=asin3φ,r≥0clear;closea=1;phy=0:pi/50:2*pi;rho=a*sin(3*phy);polar(phy,rho).三叶玫瑰线r=acos3φ,r≥0clear;closea=1;phy=0:pi/50:2*pi;rho=a*cos(3*phy);polar(phy,rho)实验二极限与导数【练习与思考】1．求下列各极限(1)lim(1-)n(2)limnn3+3n (3)lim(n+2-2n+1+n)n→∞ n→∞ n→∞clear;symsny1=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(4(4)lim(2-1)x→1x2-1x-1clear;(5)limxcot2x(6)lim(x2+3x-x)x→0 x→∞symsX;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/2y5=1/2y6=3/2m(7)lim(cos)xxm(7)lim(cos)xx→∞ xclear;(8)lim(1-1)x→1xex-1（9）limx→0xsymsxmy7=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=1y8=(exp(1)-2)∕(exp(1)-1)y9=1/3.考虑函数f(X)=3X2sin(X3),-2<X<2作出图形，并说出大致单调区间；使用diff求f)(X),并求f(X)确切的单调区间。clear;close;symsx;f=3*x^2*sin(x^3);ezplot(f,[-2,2])gridon大致的单调增区间：[-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])gridonaxis([-2.1,2.1,-60,120])f1=6*x*sin(x^3)+9*x^4*cos(x^3)用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.9948x2=-1.6926x3=-1.2401x4=0x5=1.2401x6=1.6926x7=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].对于下列函数完成下列工作，并写出总结报告，评论极值与导数的关系，(i)作出图形，观测所有的局部极大、局部极小和全局最大、全局最小值点的粗略位置；(iI)求f'(x)所有零点(即f(X)的驻点)；(iii)求出驻点处f(x)的二阶导数值；(iv)用fmin求各极值点的确切位置；(v)局部极值点与f,(x),f(x)有何关系？f(x)=X2sin(x2-X-2),X∈[-2,2]⑵f(X)=3X5-20X3+10,X∈[-3,3]f(X)=IX3-X2-X-2∣,X∈[0,3]clear;close;symsx;f=x^2*sin(x^2-x-2)ezplot(f,[-2,2])gridonf=x^2*sin(x^2-x-2)局部极大值点为：-1.6,局部极小值点为为：-0.75，-1.6全局最大值点为为：-1.6，全局最小值点为：-3f1=diff(f,x,1)ezplot(f1,[-2,2])line([-5,5],[0,0])gridonaxis([-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.5326x2=-0.7315x3=-3.2754e-027x4=1.5951ff=@(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.0272ans=2.2364ans=-0.3582ans=-9.7549e-054ans=-2.2080ans=0实验三级数【练习与思考】用taylor命令观测函数y=f(x)的Maclaurin展开式的前几项,然后在同一坐标系里作出函数y=f(x)和它的Taylor展开式的前几项构成的多项式函数的图形，观测这些多项式函数的图形向y=f(x)的图形的逼近的情况f(x)=arcsinxclear;symsxy=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=0y2=x^3/6+xy3=(35*x^9)/1152+(5*x^7)/112+(3*x^5)/40+x^3/6+xy4=

(231*x^13)∕13312+(63*x^11)∕2816(231*x^13)∕13312+(63*x^11)∕2816+(3*x^5)∕40+x^3∕6+Xf(x)=arctanxclear;symsXy=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)yl=Xy2=X-x^3∕3yff(x)=ex2clear;symsxy=exp(x^2);y1=taylor(y,0,3)y2=taylor(y,0,5)x^9∕9-x^7∕7+x^5∕5-x^3∕3+Xy4=x^13∕13-x^11∕11+x^9∕9-x^7∕7+x^5∕5-x^3∕3+X

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,':',TineWidth',3)yl=x^2+1y2=x^4∕2+x^2+1y3=x^8∕24+x^6∕6+x^4∕2+x^2+1y4=x^14∕5040+x^12∕720+x^10∕120+x^8∕24+x^6∕6+x^4∕2+x^2+1f(x)=sin2xclear;symsxy=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=0x^2-x^4∕3y2=x^2-x^4∕3y3=-x^8∕315+(2*x^6)∕45-x^4∕3+x^2y4=(4*x^14)∕42567525-(2*x^12)∕467775+(2*x^10)∕14175-x^8∕315+(2*x^6)∕45-x^4∕3+x^2f(X)=F1-Xclear;symsXy=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)yi=(5*x^2)∕2+2*x+1y2=(65*x^4)∕24+(8*x^3)∕3+(5*x^2)∕2+2*x+1y3=(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+1y4=(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

f(X)=ln(x+YI+X2)clear;symsXy=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)yl=Xy2=X-x^3∕6y3=(35*x^9)∕1152-(5*x^7)∕112+(3*x^5)∕40-x^3∕6+Xy4=(231*x^13)∕13312-(63*x^11)∕2816+(35*x^9)∕1152-(5*x^7)∕112+(3*x^5)∕40-x^3∕6+Xn2k mn=1 k.求公式∑—=吧,k=1,2,n2k mn=1 kkk=[45678];symsnsymsum(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].利用公式∑-1=e来计算e的近似值。精确到小数点后100位,这时应计算到n!n=0这个无穷级数的前多少项?请说明你的理由.解：Matlab代码为clear;clc;closeepsl=1.0e-100;ep=1;fn=1;a=1;n=1;whileep>epsla=a+fn;n=n+1;fn=fn/n;ep=fn;endfnvpa(a,100)nfn=8.3482e-101ans=2.71828182845904553488480814849026501178741455078125n=70精确到小数点后100位,这时应计算到这个无穷级数的前71项,理由是误差小于10的负100次方，需要最后一项小于10的负100次方，由上述循环知n=70时最后一项小于10的负100次方，故应计算到这个无穷级数的前71项..用练习3中所用观测法判断下列级数的敛散性∑—n2+n3n=1clear;clc;epsl=0.000001;N=50000;p=1000;symsnUn=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)ifsa<epsldisp('收敛')elsedisp('发散')end级数1∕(n^3+n^2)收敛clear;closesymsns=[];fork=1:100s(k)=symsum(1/(n^3+n^2),1,k);endplot(s,'.')⑵∑;n2nn=1clear;clc;epsl=0.000001;N=50000;p=1000;symsnUn=1/(n*2^n);SI=SymSUm(Un,1,N);s2=symsum(Un,1,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)ifsa<epsldisp('收敛')elsedisp('发散')end级数1∕(2^n*n)收敛clear;closesymsns=[];fork=1:100s(k)=symsum(1∕(2^n*n),1,k);endplot(s,'.')∑sin1nn=1clear;clc;epsl=0.00000000000001;N=50000;p=100;symsnUn=1/sin(n);SI=SymSUm(Un,1,N);s2=symsum(Un,1,N+p);sa=vpa(s2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)ifabs(sa)<epsldisp('收敛')elsedisp('发散')end级数1/sin(n)发散clear;closesymsns=[];fork=1:100s(k)=symsum(1∕sin(n),1,k);endplot(s,'.')发散∑lnnn3n=1clear;clc;epsl=0.0000001;N=50000;p=1000;symsnUn=log(n)/(n^3);SI=SymSUm(Un,1,N);s2=symsum(Un,1,N+p);Sa=VPa(S2-s1);sa=setstr(sa);sa=str2num(sa);fprintf('级数')disp(Un)ifsa<epsldisp('收敛')elsedisp('发散')end级数log(n)∕n^3收敛clear;closesymsns=[];fork=1:100s(k)=symsum(log(n)/n^3,1,k);endplot(s,'.')⑸∑史nnn=1clear;closesymsns=[];he=0;fork=1:100he=he+factorial(k)/k^k;s(k)=he;endplot(s,'.')⑹∑，(lnn)nn=3clear;clc;epsl=0.0000001;N=50000;p=1000;symsnUn=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)ifsa<epsldisp('收敛')elsedisp('发散')end级数1∕log(n)^n收敛clear;closesymsns=[];fork=3:100s(k)=symsum(1Aog(n)^n,3,k);endplot(s,'.')(7)∑nlnnn=1clear;clc;epsl=0.0000001;N=50000;p=100;symsnUn=1/(log(n)*n);SI=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)<epsldisp('收敛')elsedisp('发散')end级数1∕(n*log(n))发散clear;closesymsns=[];fork=3:300s(k)=symsum(1∕(n*log(n)),2,k);endplot(s,'.')⑻∑(-1)nnn2+1n=1clear;clc;epsl=0.0000001;N=50000;p=100;symsnUn=(-1)^n*n/(n^2+1);SI=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)<epsldisp('收敛')elsedisp('发散')end

级数((-1)^n*n)/(n^2+1)收敛clear;closesymsns=[];fork=3:300s(k)=symsum((-1)^n*n∕(n^2+1),2,k);endplot(s,'.')实验四积分【练习与思考】并用diff验证JarcsinXdx,Jsec3并用diff验证JarcsinXdx,Jsec3Xdx，JxsinX2dx,Jdx,J上,1+cosx ex+1解：Matlab代码为：symsxy1=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)/2f2=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．（定积分）2．（定积分）用trapz,quad,int计算下列定积分sinx1dx,0xJIxxdx,0J2πexsin(2X)dx,J100e-X2dx解：Matlab代码为clear;x=(0+eps):0.05:1;y1=sin(x)./x;f1=trapz(x,y1)f1=0.9460fun1=@(x)Sin(X)./x;f12=quad(fun1,0+eps,1)f12= 0.9461f13=vpa(int('sin(x)/x',0,1),5)f13=0.946083∙(椭圆的周长)用定积分的方法计算椭圆x2+止=1的周长94解：椭圆的参数方程为厂3y=2sint由参数曲线的弧长公式得S=J2π\:X'(t)2+y'(t)2dt=J2πJ9sin21+4cos2tdt=J2π√5sin21+4dt00 0Matlab代码为s=vpa(int('sqrt(5*sin(t)^2+4)','t',0,2*pi),5)s=15.865.(二重积分)计算数值积分∫∫(1+X+y)dxdyX2+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或符号积分法：symsXyxi=int(1+x+y,y,1-sqrt(1-x^2),1+sqrt(1-x^2));s=int(xi,x,-1,1)s=2*pi.(假奇异积分)用trapz,quad8计算积分∫1x1/3cosXdX,会出现什么问题？-1分析原因，并求出正确的解。解：Matlab代码为clearx=-1:0.05:1;y=x.^(1/3).*cos(x);s1=trapz(x,y)fun5=@(x)xC(1/3).*cos(x);s2=quad(fun5,-1,1)int('x^(1/3)*cos(x)','x',-1,1)s1=0.9036+0.5217is2=0.9114+0.5262iWarning:Explicitintegralcouldnotbefound.ans=int(x^(1/3)*cos(x),x=-1..1)，原函数不存在，不能用int函数运算。用梯形法和辛普森法计算数值积分时，由于对负数的开三次方运算结果为复数，所以导致结果错误且为复数；显然被积函数为奇函数，在对称区间上的积分等于0，此时可以这样处理：(1)重新定义被积函数%fun5.mfunctiony=fun5(x)[m,n]=size(x);fork=1:mforl=1:ny(k,l)=nthroot(x(k,l),3)*cos(x(k,l));endendend用辛普森法：s=quad('fun5',-1,1)s=0用梯形法clear;x=-1:0.01:1;y=fun5(x);s=trapz(x,y)s=-1.3878e-0176∙(假收敛现象)考虑积分I(k)=∫kπ∣sinXldx,0(1)用解析法求I(k);clear;symsxk;Ik=int(abs(sin(x)),0,k*pi)Warning:Explicitintegralcouldnotbefound.Ik=int(abs(sin(x)),x=0..pi*k)(2)分别用trapz,quad和quad8求I(4),I(6)和I(8),发现什么问题？clear;fork=4:2:8;x=0:pi/1000:k*pi;y=abs(sin(x));trapz(x,y)endans=8.0000ans=12.0000ans=16.0000fork=4:2:8fun6=@(x)abs(sin(x));quad(fun6,0,k*pi)endans=8.0000ans=12.0000ans=16.000016.00007．(Simpson积分法)编制一个定步长Simpson法数值积分程序.计算公式为h. -… - … - 、I≈S=-(f+4f+2f+4f+…+2f+4f+f)

n31 2 3 4 n-1 n n+1b-a其中n为偶数,-= ,f=f(a+(i-1)-),i=1,2,…,n+1.ni解：Matlab代码为%fun7.mfunctiony=fun7(f_name,a,b,n)%f_name为被积函数%[a,b]为积分区间%n为偶数，用来确定步长h=(b-a)/nifmod(n,2)~=0disp('n必须为偶数')return;endifnargin<4n=100;endifnargin<3disp('请输入积分区间')endifnargin==0disp('error')endh=(b-a)/n;x=a:h:b;s=0;fork=1:n+1ifk==1||k==(n+1)xishu=1;elseifmod(k,2)==0xishu=4;elsexishu=2;ends=s+feval(f_name,x(k))*xishu;endy=s*h/3;end8．(广义积分)计算广义积分exp(-X2) tan(x) SinXdx,J——^―dx,J dx-∞1+X4 0√X 0√1-X2并验证公式exp(- )∫∞ 2dx=1,∫∞Sinxdx=∏.-∞ V2π 0X 2解：Matlab代码为clear;symsXSI=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.4348s2=0.7968s3=0.8933s4=1.0s5=pi/2-sinint(1/4503599627370496)实验五二元函数的图形【练习与思考】.画出空间曲线Z=IOSin、；K在-30<X,y<30范围内的图形，并画v'1+X2+y2出相应的等高线。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(I+X∙^2+Y.^2);mesh(X,Y,Z)Sia)Sia).根据给定的参数方程，绘制下列曲面的图形。椭球面x=3cosusinv,y=2cosucosv,z=sinuclear;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=4u2clear;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)axisequal3secusinv,y2secucosv,z4tanuclear;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)双曲抛物面X=Ud)双曲抛物面X=U,y=V,Z=U2-V23clearu=-3:0.1:3;[U,V]=meshgrid(u);x=U;y=V;z=(U.^2-V.^2)∕3;mesh(x,y,z)e)旋转面x=e)旋转面x=lnusinv,y=lnucosv,z=uclear;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)axisequalf)圆锥面xf)圆锥面x=usinv,y=ucosv,z=uclear;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)axisequalZ=0.4SinVu=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)yucosv,z4vclear;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,试拟合一曲面,确定合适的模型,并由此找出最高点和该点的高程.表5-2高程数据yx100200300400100636697624478200698712630478300680674598412400662626552334clc;clear;x1=[100100100100200200200200300300300300400400400400];x2=[100200300400100200300400100200300400100200300400];y=[636698680662697712674626624630598552478478412334]';x=[x1',x2'];x0=[11111];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)Localminimumpossible.lsqcurvefitstoppedbecausethefinalchangeinthesumofsquaresrelativetoitsinitialvalueislessthanthedefaultvalueofthefunctiontolerance.

beta=Columns1through3538.4375 1.4901 0.6189Columns4through5-0.0046 -0.0017contour(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:OptionsLargeScale='off'andAlgorithm='trust-region-reflective'conflict.IgnoringAlgorithmandrunningactive-setalgorithm.Toruntrustregion-reflective,setLargeScale='on'.Torunactive-setwithoutthiswarning,useAlgorithm='active-set'.>Infminconat445Localminimumpossible.Constraintssatisfied.fminconstoppedbecausethepredictedchangeintheobjectivefunctionislessthanthedefaultvalueofthefunctiontoleranceandconstraintsweresatisfiedtowithinthedefaultvalueoftheconstrainttolerance.Noactiveinequalities.x=161.9676182.0320fval=-715.4403

heigh和height两个函数分别定义如下：(应写在m文件中)%heigh.mfunctionf=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.mfunctiony=height(x)y=-(538.4375+1.4901*x(1)+0.6189*x(2)-0.0046*x(1).^2-0.0017*x(2).^2);end实验六多元函数的极值【练习与思考】.求Z=X4+y4-4与+1的极值，并对图形进行观测。解：Maltab代码为symsXy;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)不是极值点。clearsymsXy;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^2B=-4C=12*y^2由判别法可知(1,1)、(-1,-1)均为极小值点。.求函数fQ,y)=X2+2y2在圆周X2+y2=1的最大值和最小值。解：构造Lagrange函数L(X,y)=X2+2y2+λ(X2+y2-1)求Lagrange函数的自由极值.先求L关于X,y,λ的一阶偏导数，再解正规方程可得所求的极值点，Matlab代码为clear;symsxykL=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)为函数的条件极值点，经判断函数fQ,y)=X2+2y2在(1,0)、(-1,0)取得极小值，在(0,1)、(0,-1)取得极大值。.在球面X2+w+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)满足X2+产+Z2=1；构造Lagrange函数L(X,y,z,λ)=(X-3)2+(y-1)2+(Z+1)2+λ(X2+产+Z2-1)求Lagrange函数的自由极值.先求L关于X,y,Z,λ的一阶偏导数，再解正规方程可得所求的极值点，Matlab代码为clearclear;symsxyzkL=(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]得到条件极值点为[5.3668,18.633]得到条件极值点为（3111,中,-斗）、（-九11,-卫1,色1）,经判断，球面X2+W+Z2=1上与点（3,1,-1）距离最近的点为（午1,]!1,-：1¥），最远的点.求函数f(X,y,Z)=X+2y+3Z在平面X-y+Z=1与柱面X2+W=1的交线上的最大值。解：构造Lagrange函数L(X,y,z,λ)=X+2y+3z+λ(X-y+Z-1)+λ(X2+y2-1)求Lagrange函数的自由极值.先求L关于XJ,Z,λ,λ的一阶偏导数，再解正规方程12可得所求的极值点，Matlab代码为clear;symsXyZklk2L=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=SoIVe(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]经判断可知，函数f(X,y,Z)=X+2y+3Z在平面X-y+Z=1与柱面X2+y2=1的交线上的最大值为3+、/29。.求函数Z=X2+产在三条直线X=1,y=1,X+y=1所围区域上的最大值和最小值。解：显然此函数的驻点为0,0)不在此区域内，因此该函数的最大值和最小值点应在三条边界上，下面分别求此函数在这三条边界上的最大值和最小值atlab代码如下(1)求函数在直线边界c=1,0≤y≤1上的最大值和最小值将x=1代入原函数，则二元函数变为一元函数z=1+y2,0≤y≤1最大值点为y=1，最大值为2,最小值点为y=0,最小值为1；(2)求函数在直线边界z=1,0≤x≤1上的最大值和最小值将x=1代入原函数，则二元函数变为一元函数z=1+x2,0≤x≤1最大值点为x=1，最大值为2,最小值点为＜=0,最小值为1；(3)求函数在直线边界＜+y=1,0≤x≤1上的最大值和最小值将y=1-x代入原函数，则二元函数变为一元函数z=(1-x)2+x2,0≤x≤1用Matla喻令求此函数的最大和最小值点先求驻点clear;symsxz=(1-x)^2+x^2;dzx=diff(z,x);x=solve(dzx,x)＜=1∕2z1=eval(z)计算在驻点处的函数值z1=1∕2计算在区间端点处的函数值z2=subs(z,0)z3=subs(z,1)z2= 1z3= 1比较函数在各点处的函数值可知函数的最大值点为(1,1)，对应的最大值为2，最小值点为(1/2,1/2)，最小值为1/2。实验七常微分方程【练习与思考】1．求下列微分方程的解析解a)一阶线性方程y，-x3y=2dsolve('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'-χy2-y=0dsolve('Dy-x*y^2-y=0','x')ans=0exρ(x)∕(C4-exρ(x)*(x-1))C)高阶线性齐次方程y'-y-3y+2y=0dsolve('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-3卢2y=3sinXdsolve('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)欧拉方程X3y'+X2y-3Xy=4X2dsolve('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=2XHy(0)=1,歹(0)=3的解析解和数值解,并进行比较解：解析解为y=dsolve('(1+x^2)*D2y=2*x*Dy','y(0)=1','Dy(0)=3','x')yx*(x^2+3)+1数值解：设y=y,y=y'则原方程化为微分方程组12•y'=y2y=—21+X2定义函数m文件fun7_2.m如下functionf=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')holdonezplot('x*(x^2+3)+1',[0,5])3．分别用ode45和ode15s求解Van-del-Pol方程d2x 、dx—-ιooqι-X2)—-X=0Vdt2 dtx(0)=0,X'(0)=1)的数值解,并进行比较.解：设X=X,X=X'则原方程化为微分方程组12X'=X12X'=1000(1-X2)X+X2 12 1定义函数m文件fun7_3.m如下functionf=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．(单摆运动的近似解析解)当单摆初始角度θ较小时,θ(≤θ)也较小,从00sinθ≈θ,单摆运动微分方程可近似写为mlθ=mgθ,θ(0)=θ,θ'(0)=00求此方程的解析解,并与练习3中的数值解进行比较.解：用Matlab命令求此方程的解析解，按练习3中的取值l=1,g=9.8,θ(0)=15clear;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.mfunctionf=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')holdonplot(t,s)

xlabel('t'),ylabel('y2')显然，在这个区间内，二者差别较大，下面改为较小区间clear;close;s=dsolve('D2y=9.8*y','y(0)=1*pi/180','Dy(O)=0','t')[t,y]=ode45('fun7_4',[0,2],[1*pi/180,0]);s=eval(s);plot(t,y(:,1),'ro');holdonplot(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］内二者能够较好吻合。实验八平面图形的几何变换【练习与思考】.将函数y=eT2的图形向右平移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,':',TineWidth',3);axis(［-3,6,-1,5］)xlabel('x');ylabel('y');gridon

.将函数y=eT2的图形在水平方向收缩一倍，在垂直方向放大一倍。clear;close;x=-2:0.1:2;y=exp(-x.^2);χ1=χ∕2;%图形在水平方向收缩一倍；y1=y*2；%图形在垂直方向放大一倍plot(x,y,x1,y,'-.',x1,y1,':',TineWidth',3);axis([-3,3,-1,3])xlabel('x');ylabel('y');gridon.将函数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:',TineWidth',3);Iegend('原图；顺时针旋转30度角后的图')xlabel('x');ylabel('y');gridon.已知函数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');gridon.做怎样的变换才能使函数图形绕给定的点Qb)转动?这个变换可以分解成3个基本变换:平移量为(-a,-b)的平移变换T,旋转角度为a的旋转变换T,T的逆变换T.1.求出变换矩阵,写出与变换相应的方程,并对具体的函数图形进行变换.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=[10-a;01-b;001];T2=[cos(alpha)-sin(alpha)0;sin(alpha)cos(alpha)0;001];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');gridon

clear;close;a=1;b=2;alpha=60*pi∕180;T1=[10-a;01-b;001];T2=[cos(alpha)-sin(alpha)0;sin(alpha)cos(alpha)0;001];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');gridon实验九π的近似计算【练习与思考】.利用勾股定理推导在割圆术中给出的公式X=、,：2-；4-X2,S=3∙2nx,S<2S-S。6-2n+1 tY6-2n6-2n+! 6-2n 2n2nn解：略.利用韦达公式，构造出一种迭代算法来计算π的近似值，并进行实际计算，评价算法效果。解：韦达(VieTa)公式TOC\o"1-5"\h\z2_√2v2+72\；2+v2+√2v2+%2+√2+、2 •π2 2 2 22s设a=v2,a=.,2+a,(n=2,3,4,…)，S=2,S=—n-ι,(n=1,2,3∙,∙∙),则lιmS=πn、 n-1 0 n nn→∞n于是得到一种迭代算法，实际计算的Matlab代码为clear;a=sqrt(2);s=2;n=15;fork=1:ns=2*s∕a;a=sqrt(2+a);pai=vpa(sym(s),20)error=vpa(sym(pi)-sym(s),20)endpai=2.8284271247461900976error=0.31316552884360314086pai=3.0614674589207178101error=0.080125194669075428318pai=3.1214451522580519693error=0.020147501331741269122pai=3.136548490545938872error=0.005044163043854366431pai=3.1403311569547525117error=0.0012614966350407267397pai=3.1412772509327724357error=0.00031540265702080273067pai=3.1415138011443008992error=0.000078852445492339280772pai=3.1415729403670913378error=0.000019713222701900686015pai=3.1415877252771595707error=0.000004928312633667782247p