數(shù)學(xué)實驗課后習(xí)題解答

數(shù)學(xué)實驗課后習(xí)題解答王汝軍編寫實驗一曲線繪圖【練習(xí)與思考】畫出下列常見曲線的圖形。以直角坐標(biāo)方程表示的曲線：.立方曲線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以參數(shù)方程表示的曲線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

以極坐標(biāo)方程表示的曲線：.阿基米德線r=aφ,r≥0clear;a=1;phy=0:pi/50:6*pi;rho=a*phy;polar(phy,rho,'r-*').對數(shù)螺線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)實驗二極限與導(dǎo)數(shù)【練習(xí)與思考】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.考慮函數(shù)f(X)=3X2sin(X3),-2<X<2作出圖形，并說出大致單調(diào)區(qū)間；使用diff求f)(X),并求f(X)確切的單調(diào)區(qū)間。clear;close;symsx;f=3*x^2*sin(x^3);ezplot(f,[-2,2])gridon大致的單調(diào)增區(qū)間：[-2,-1.7],[-1.3,1.2],[1.7,2];大致的單點減區(qū)間：[-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函數(shù)找f'(X)的零點，即原函數(shù)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確切的單調(diào)增區(qū)間:[-1.9948,-1.6926],[-1.2401,1.2401],[1.6926,1.9948]確切的單調(diào)減區(qū)間:[-2,-1.9948],[-1.6926,-1.2401],[1.2401,1.6926],[1.9948,2].對于下列函數(shù)完成下列工作，并寫出總結(jié)報告，評論極值與導(dǎo)數(shù)的關(guān)系，(i)作出圖形，觀測所有的局部極大、局部極小和全局最大、全局最小值點的粗略位置；(iI)求f'(x)所有零點(即f(X)的駐點)；(iii)求出駐點處f(x)的二階導(dǎo)數(shù)值；(iv)用fmin求各極值點的確切位置；(v)局部極值點與f,(x),f(x)有何關(guān)系？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函數(shù)找f'(x)的零點，即原函數(shù)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實驗三級數(shù)【練習(xí)與思考】用taylor命令觀測函數(shù)y=f(x)的Maclaurin展開式的前幾項,然后在同一坐標(biāo)系里作出函數(shù)y=f(x)和它的Taylor展開式的前幾項構(gòu)成的多項式函數(shù)的圖形，觀測這些多項式函數(shù)的圖形向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的近似值。精確到小數(shù)點后100位,這時應(yīng)計算到n!n=0這個無窮級數(shù)的前多少項?請說明你的理由.解：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精確到小數(shù)點后100位,這時應(yīng)計算到這個無窮級數(shù)的前71項,理由是誤差小于10的負100次方，需要最后一項小于10的負100次方，由上述循環(huán)知n=70時最后一項小于10的負100次方，故應(yīng)計算到這個無窮級數(shù)的前71項..用練習(xí)3中所用觀測法判斷下列級數(shù)的斂散性∑—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('級數(shù)')disp(Un)ifsa<epsldisp('收斂')elsedisp('發(fā)散')end級數(shù)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('級數(shù)')disp(Un)ifsa<epsldisp('收斂')elsedisp('發(fā)散')end級數(shù)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('級數(shù)')disp(Un)ifabs(sa)<epsldisp('收斂')elsedisp('發(fā)散')end級數(shù)1/sin(n)發(fā)散clear;closesymsns=[];fork=1:100s(k)=symsum(1∕sin(n),1,k);endplot(s,'.')發(fā)散∑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('級數(shù)')disp(Un)ifsa<epsldisp('收斂')elsedisp('發(fā)散')end級數(shù)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('級數(shù)')disp(Un)ifsa<epsldisp('收斂')elsedisp('發(fā)散')end級數(shù)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('級數(shù)')disp(Un)if(sa)<epsldisp('收斂')elsedisp('發(fā)散')end級數(shù)1∕(n*log(n))發(fā)散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('級數(shù)')disp(Un)if(sa)<epsldisp('收斂')elsedisp('發(fā)散')end

級數(shù)((-1)^n*n)/(n^2+1)收斂clear;closesymsns=[];fork=3:300s(k)=symsum((-1)^n*n∕(n^2+1),2,k);endplot(s,'.')實驗四積分【練習(xí)與思考】并用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解：橢圓的參數(shù)方程為廠3y=2sint由參數(shù)曲線的弧長公式得S=J2π\(zhòng):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.(二重積分)計算數(shù)值積分∫∫(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,會出現(xiàn)什么問題？-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)，原函數(shù)不存在，不能用int函數(shù)運算。用梯形法和辛普森法計算數(shù)值積分時，由于對負數(shù)的開三次方運算結(jié)果為復(fù)數(shù)，所以導(dǎo)致結(jié)果錯誤且為復(fù)數(shù)；顯然被積函數(shù)為奇函數(shù)，在對稱區(qū)間上的積分等于0，此時可以這樣處理：(1)重新定義被積函數(shù)%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?(假收斂現(xiàn)象)考慮積分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),發(fā)現(xiàn)什么問題？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法數(shù)值積分程序.計算公式為h. -… - … - 、I≈S=-(f+4f+2f+4f+…+2f+4f+f)

n31 2 3 4 n-1 n n+1b-a其中n為偶數(shù),-= ,f=f(a+(i-1)-),i=1,2,…,n+1.ni解：Matlab代碼為%fun7.mfunctiony=fun7(f_name,a,b,n)%f_name為被積函數(shù)%[a,b]為積分區(qū)間%n為偶數(shù)，用來確定步長h=(b-a)/nifmod(n,2)~=0disp('n必須為偶數(shù)')return;endifnargin<4n=100;endifnargin<3disp('請輸入積分區(qū)間')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)實驗五二元函數(shù)的圖形【練習(xí)與思考】.畫出空間曲線Z=IOSin、；K在-30<X,y<30范圍內(nèi)的圖形，并畫v'1+X2+y2出相應(yīng)的等高線。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).根據(jù)給定的參數(shù)方程，繪制下列曲面的圖形。橢球面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)旋轉(zhuǎn)面x=e)旋轉(zhuǎn)面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高程數(shù)據(jù)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');%設(shè)置下界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兩個函數(shù)分別定義如下：(應(yīng)寫在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實驗六多元函數(shù)的極值【練習(xí)與思考】.求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]經(jīng)計算可知，函數(shù)的駐點為(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)均為極小值點。.求函數(shù)fQ,y)=X2+2y2在圓周X2+y2=1的最大值和最小值。解：構(gòu)造Lagrange函數(shù)L(X,y)=X2+2y2+λ(X2+y2-1)求Lagrange函數(shù)的自由極值.先求L關(guān)于X,y,λ的一階偏導(dǎo)數(shù)，再解正規(guī)方程可得所求的極值點，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)為函數(shù)的條件極值點，經(jīng)判斷函數(shù)fQ,y)=X2+2y2在(1,0)、(-1,0)取得極小值，在(0,1)、(0,-1)取得極大值。.在球面X2+w+Z2=1求出與點(3,1,-1)距離最近和最遠點。解：設(shè)球面上的點為(x,y,z),則此點與點(3,1,-1)的距離為d(X,y,Z)=(X-3)2+(y-1)2+(Z+1)2且(x,y,z)滿足X2+產(chǎn)+Z2=1；構(gòu)造Lagrange函數(shù)L(X,y,z,λ)=(X-3)2+(y-1)2+(Z+1)2+λ(X2+產(chǎn)+Z2-1)求Lagrange函數(shù)的自由極值.先求L關(guān)于X,y,Z,λ的一階偏導(dǎo)數(shù)，再解正規(guī)方程可得所求的極值點，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,-衛(wèi)1,色1）,經(jīng)判斷，球面X2+W+Z2=1上與點（3,1,-1）距離最近的點為（午1,]!1,-：1￥），最遠的點.求函數(shù)f(X,y,Z)=X+2y+3Z在平面X-y+Z=1與柱面X2+W=1的交線上的最大值。解：構(gòu)造Lagrange函數(shù)L(X,y,z,λ)=X+2y+3z+λ(X-y+Z-1)+λ(X2+y2-1)求Lagrange函數(shù)的自由極值.先求L關(guān)于XJ,Z,λ,λ的一階偏導(dǎo)數(shù)，再解正規(guī)方程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]經(jīng)判斷可知，函數(shù)f(X,y,Z)=X+2y+3Z在平面X-y+Z=1與柱面X2+y2=1的交線上的最大值為3+、/29。.求函數(shù)Z=X2+產(chǎn)在三條直線X=1,y=1,X+y=1所圍區(qū)域上的最大值和最小值。解：顯然此函數(shù)的駐點為0,0)不在此區(qū)域內(nèi)，因此該函數(shù)的最大值和最小值點應(yīng)在三條邊界上，下面分別求此函數(shù)在這三條邊界上的最大值和最小值atlab代碼如下(1)求函數(shù)在直線邊界c=1,0≤y≤1上的最大值和最小值將x=1代入原函數(shù)，則二元函數(shù)變?yōu)橐辉瘮?shù)z=1+y2,0≤y≤1最大值點為y=1，最大值為2,最小值點為y=0,最小值為1；(2)求函數(shù)在直線邊界z=1,0≤x≤1上的最大值和最小值將x=1代入原函數(shù)，則二元函數(shù)變?yōu)橐辉瘮?shù)z=1+x2,0≤x≤1最大值點為x=1，最大值為2,最小值點為＜=0,最小值為1；(3)求函數(shù)在直線邊界＜+y=1,0≤x≤1上的最大值和最小值將y=1-x代入原函數(shù)，則二元函數(shù)變?yōu)橐辉瘮?shù)z=(1-x)2+x2,0≤x≤1用Matla喻令求此函數(shù)的最大和最小值點先求駐點clear;symsxz=(1-x)^2+x^2;dzx=diff(z,x);x=solve(dzx,x)＜=1∕2z1=eval(z)計算在駐點處的函數(shù)值z1=1∕2計算在區(qū)間端點處的函數(shù)值z2=subs(z,0)z3=subs(z,1)z2= 1z3= 1比較函數(shù)在各點處的函數(shù)值可知函數(shù)的最大值點為(1,1)，對應(yīng)的最大值為2，最小值點為(1/2,1/2)，最小值為1/2。實驗七常微分方程【練習(xí)與思考】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的解析解和數(shù)值解,并進行比較解：解析解為y=dsolve('(1+x^2)*D2y=2*x*Dy','y(0)=1','Dy(0)=3','x')yx*(x^2+3)+1數(shù)值解：設(shè)y=y,y=y'則原方程化為微分方程組12?y'=y2y=—21+X2定義函數(shù)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)的數(shù)值解,并進行比較.解：設(shè)X=X,X=X'則原方程化為微分方程組12X'=X12X'=1000(1-X2)X+X2 12 1定義函數(shù)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．(單擺運動的近似解析解)當(dāng)單擺初始角度θ較小時,θ(≤θ)也較小,從00sinθ≈θ,單擺運動微分方程可近似寫為mlθ=mgθ,θ(0)=θ,θ'(0)=00求此方程的解析解,并與練習(xí)3中的數(shù)值解進行比較.解：用Matlab命令求此方程的解析解，按練習(xí)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練習(xí)3的中數(shù)值解%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')顯然，在這個區(qū)間內(nèi)，二者差別較大，下面改為較小區(qū)間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

由圖象可知，區(qū)間改為［0,2］上時能觀察出大概在區(qū)間［0,1.5］內(nèi)二者能夠較好吻合。實驗八平面圖形的幾何變換【練習(xí)與思考】.將函數(shù)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

.將函數(shù)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.將函數(shù)y=X2的圖形以原點為中心，順時針旋轉(zhuǎn)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('原圖；順時針旋轉(zhuǎn)30度角后的圖')xlabel('x');ylabel('y');gridon.已知函數(shù)y=2x-x2,,0≤x≤2，試擴展函數(shù)的定義域，使之成為以2周期的偶函數(shù)，并畫出函數(shù)在［-8，8］上的圖形。若要把函數(shù)延拓成以4為周期的奇函數(shù)呢?解：延拓成以2周期的偶函數(shù)，畫出函數(shù)在［-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');把函數(shù)延拓成以4為周期的奇函數(shù)，畫出函數(shù)在［-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.做怎樣的變換才能使函數(shù)圖形繞給定的點Qb)轉(zhuǎn)動?這個變換可以分解成3個基本變換:平移量為(-a,-b)的平移變換T,旋轉(zhuǎn)角度為a的旋轉(zhuǎn)變換T,T的逆變換T.1.求出變換矩陣,寫出與變換相應(yīng)的方程,并對具體的函數(shù)圖形進行變換.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í)與思考】.利用勾股定理推導(dǎo)在割圓術(shù)中給出的公式X=、,：2-；4-X2,S=3?2nx,S<2S-S。6-2n+1 tY6-2n6-2n+! 6-2n 2n2nn解：略.利用韋達公式，構(gòu)造出一種迭代算法來計算π的近似值，并進行實際計算，評價算法效果。解：韋達(VieTa)公式TOC\o"1-5"\h\z2_√2v2+72\；2+v2+√2v2+%2+√2+、2 ?π2 2 2 22s設(shè)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