Matlab 实验程序 汇总40

MATLAB数学实验程序(2011-12-18)ex1.单位圆的正多边形逼近alpha=linspace(0,2*pi,13);x=sin(alpha);y=cos(alpha);I2=1:2:13;I1=1:4:13;x1=x(I1);y1=y(I1);x2=x(I2);y2=y(I2);plot(x1,y1,x2,y2,x,y)axisoffex2抛射曲线的最大射程计算v0=100;g=9.8;alpha=pi/4;T=2*v0*sin(alpha)/g;t=linspace(0,T,17);x=v0*t*cos(alpha);y=vO*t*sin(alpha)-g*t42/2;plot(x,y,x,y,'r*')xmax=x(end)ex3.巴拿马草帽.图形[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.A2+y.A2)+eps;z=sin(r)./r;mesh(x,y,z)colormap([1,0,0])axisoffex4.三叶玫瑰线实验theta=linspace(0,2*pi,10000);n=3;r=cos(n*theta);figure(1)polar(theta,r,'k')x=r.*cos(theta);y=r.*sin(theta);figure(2)comet(x,y)ex5.抛射曲线簇绘图程序alpha=linspace(0.1,pi/2,20);v0=100;g=9.8;Taim=2*v0*sin(alpha)/g;T=(0:16)'*Taim/16;X=v0*T*diag(cos(alpha));Y=v0*T*diag(sin(alpha))-g*T.A2/2;%计算正十二边形顶点坐标%提取正三边形坐标%提取正六边形坐标%设定初速度%计算飞行时间%设定飞行时刻%计算坐标点%显示最远射程%创建网格点矩阵%计算网格点到原点距离%计算二元函数值%二元函数绘图%设置红色%去掉坐标架%创建区间［0，2］上离散点%设置三叶玫瑰线参数%计算玫瑰线坐标数据%创建第一个图形窗%极坐标绘图%直角坐标转换%创建第二个图形窗口%用慧星绘图做动态模拟%创建20个发射角数据的行向量%计算对应的飞行时间%计算飞行时刻矩阵T%对角矩阵右乘飞行时刻矩阵%%创建网格点%设置复数自变量%计算复变函数值%分离实部和虚部plot(X,Y,'r')ex6.黎曼曲面绘图程序[r,theta]=meshgrid(s,t);z=r.*exp(i*theta);f=sqrt(r).*exp(i*theta/2);x=real(z);y=imag(z);figure(1),mesh(x,y,real(f)),axisofffigure(2),mesh(x,y,imag(f)),axisoffpk=P(ki,:);L0=a;forkj=1:3L=L0;II=pk(kj);L(II)=[];pk=P(ki,:);L0=a;forkj=1:3L=L0;II=pk(kj);L(II)=[];A(ki,kj)=prod(L);L0(II)=b(II);endcolormap([001]),view(AZ,20)M(:,k)=getframe;AZ=AZ+15;endt=linspace(-pi/4,pi/4,50);ro=sqrt(cos(2*t));x=ro.*cos(t);y=ro.*sin(t);X=[0.6;1]*x;Y=[0.6;1]*y;Z=[0;1]*ones(1,50);mesh(X,Y,Z),holdoncolormap([000])view(17,-12),axisofffill3(0.6*x,0.6*y,zeros(1,50),'c')exp7.宝石切割问题实验程序P=perms(l:3);a=[19,14,12];b=[5,4,2];forki=1:6endQ=A*[1;1;1];disp([P,2*Q])exp8马鞍面动态旋转帧播放程序M=moviein(16);[x,y]=meshgrid(-6:.5:6);z=x.A2-y.A2;mesh(x,y,z),axisoffAZ=-30;fork=1:24movie(M,2)ex9.飞机机翼图仿真程序%创建全排歹【I操作数%设置石料数据%提取操作数%模拟切割%计算单次切割面积%计算每种操作切割面积%创建马鞍面坐标数据%马鞍面绘图

%设置初始方位角%旋转图形%帧设计%方位解增加15。%帧播放%设置50个点%计算双纽线坐标%直角坐标转换%XY坐标扩充矩阵%Z坐标矩阵%设置黑色%设置视角、去坐标架%底面填充Y(k,2:k-1)=Y(k-1,1:k-2)+Y(k-1,2:k-1);endY(k,2:k-1)=Y(k-1,1:k-2)+Y(k-1,2:k-1);endex13.正方形旋转实验程序%输入正整数%n不等于l%求n/2的余数%第一种操作%第二种操作exlO.角谷猜想或3n+l问题的周期计算n=input('inputn:=');N=n;T=l;whilen~=lr=rem(n,2);ifr==0

n=n/2;elsen=3*n+l;endN=[N,n];T=T+l;endT,NExll.骑摩托问题模拟程序A=0;B=l00;va=l0;vb=8;vc=60;f=l;k=0;while(B-A)>0.2iff==ltk=(B-A)/(vb+vc);elsetk=(B-A)/(vc+va);endA=A+va*tk;B=B-vb*tk;f=-f;k=k+1;endk,A,Bex12.杨辉三角形计算程序functionY=young(n)ifnargin==0,n=3;endY=eye(n);Y(:,1)=ones(n,1);fork=3:n%计算C和B相遇时间%计算C和A相遇时间%根据相遇时间计算位移%将矩阵第一列和对角元置为1%根据前一行计算当前行数据xy=[-1-1;1-1;11;-11;-1-1];A=[cos(pi/24)—sin(pi/24);sin(pi/24)cos(pi/24)];%创建正交矩阵x=xy(:,1);y=xy(:,2); %提取坐标数据%画线并暂停一秒%旋转并压缩%画线并暂停一秒%旋转并压缩fork=1:24xy=.89*xy*A';x=xy(:,1);y=xy(:,2);line(x,y),pause(1)endex14.中国农历年与十二生肖functionyears(year)%输入年份year%输出天干地支农历年生肖属相

ifnargin==0,year=2011;endS1='甲乙丙丁戊己庚辛壬癸ifnargin==0,year=2011;endS1='甲乙丙丁戊己庚辛壬癸'；S2='子丑寅卯辰巳午未申酉戍亥';S='鼠牛虎兔龙蛇马羊猴鸡狗猪'；k1=mod(year-4,10)+1；k2=mod(year-4,12)+1;s1=S1(k1);s2=S2(k2);s=S(k2);X=[S1,'..';S2;S]disp([int2str(year)，年是',sl,s2,'年,%创建天干字符串数组%创建地支字符串数组%创建十二生肖字符串数组%年份除10的余数%准确定位%输出天干地支十二生肖表','属相：’,s])ex15.生日与十二星座模拟functionAries(x)%输入生日格式[mm,dd]data=[1,19;2,18;3,20;4,20;5,20;6,21;7,23;8,23;9,23;10,24;11,22;12,22];Const=['水瓶';'双鱼';'白羊';'金牛';'双子';'巨蟹';'狮子';'处女';'天秤';'天蝎';'射手';'魔蝎'];P=data(:,1);Q=data(:,2);p=x(1);q=x(2);II=find(P<=p&p<P+1);ifq<Q(II),II=II-1;endifII==0,II=12;endzodiac=Const'disp(['你输入的生日属于:',Const(II,:),座'])Ex16科赫曲线绘图程序n=1;p=[00;100];A=[cos(pi/3)-sin(pi/3);sin(pi/3)cos(pi/3)]; %创建正交矩阵fork=1:5j=1;fori=1:nq1=p(i,:);q2=p(i+1,:);fori=1:nq1=p(i,:);q2=p(i+1,:);d=(q2-q1)/3;r(j,:)=q1;r(j+1,:)=q1+d;r(j+2,:)=q1+d+d*A';r(j+3,:)=q1+2*d;j=j+4;%提取线段起点%提取线段终点%计算第一个插入点%计算第二个插入点%计算第三个插入点endn=4*n;p=[];p=[r;q2];%整合分形曲线的数据endx=p(:,1);y=p(:,2);plot(x,y)%分离出坐标数据ex17•endn=4*n;p=[];p=[r;q2];%整合分形曲线的数据endx=p(:,1);y=p(:,2);plot(x,y)%分离出坐标数据ex17•北京到纽约飞行航程计算程序city=[40,116;31,122;36,140;37,-123;41,-76];R=6400+10;theta=city(:,1)*pi/180;fai=city(:,2)*pi/180;x=R*cos(theta).*cos(fai);y=R*cos(theta).*sin(fai);%五城市经纬度%飞行高度设置%提取纬度数据%提取经度数据%球坐标变换%城市坐标%城市坐标%计算距离矩阵op=[x,y,z];formatbankDmatrix=R*acos(op*op7R^2)Ex18.北京到纽约飞行模拟实验程序functionflytravel()city=[40,116;31,122;36,140;37,-123;41,-76];figure(1),sphere(24),colormap([111])axisoff,holdonp1=city(1,:);p2=city(5,:);skyway(p1,p2)figure(2),sphere(24),colormap([1,1,1])axisoff,holdonfork=1:4p1=city(k,:);p2=city(k+1,:);skyway(p1,p2)endfunctionskyway(p1,p2)city=[p1;p2]*pi/180;theta=city(:,1);fai=city(:,2);x=cos(theta).*cos(fai);y=cos(theta).*sin(fai);z=sin(theta);t=linspace(0,1,20);xt=(1-t)*x(1)+t*x(2);yt=(1-t)*y(1)+t*y(2);zt=(1-t)*z(1)+t*z(2);r=sqrt(xt.*xt+yt.*yt+zt.*zt);xt=xt./r;yt=yt./r;zt=zt./r;plot3(x,y,z,'ro',xt,yt,zt,'b','linewidth',2)exl9•牟合方盖体积计算与绘图程序symsxyRy0=sqrt(RA2-xA2);V=8*int(int('sqrt((RA2-xA2))',y,0,y0),x,0,R)V0=subs(V,R,l)r=linspace(0,l,20);t=linspace(0,2*pi,80);u=r'*cos(t);v=r'*sin(t);w=sqrt(1+eps-u.A2);meshz(u,v,w),colormap([0,0,l])axisoffview(-46,52)colormap([0,0,0])ex20•左矩形公式计算及绘图程序

functionS=interg(h)ifnargin==0,h=.5;endf=inline('x.A3./(exp(x)-l)');x=eps:h:5-h;y=f(x);n=length(x);X(l:3:3*n)=x;X(2:3:3*n)=x;X(3:3:3*n)=[x(2:n),5];Y(2:3:3*n)=y;Y(3:3:3*n)=y;X=[X,5];Y=[Y,0];t=eps:0.l:l0;u=f(t);plot(X,Y,'k',t,u,'k','LineWidth',l)S=sum(y)*h;Ex2l.探月卫星初始轨道的最大速度计算程序h=200;H=5l000;R=6378;a=(h+H+2*R)/2;c=(H-h)/2;el=c/a;b=sqrt(a*a-c*c);T=l6*3600;Sl=a*b*pi/T;symse2th=200;H=5l000;R=6378;a=(h+H+2*R)/2;c=(H-h)/2;el=c/a;b=sqrt(a*a-c*c);T=l6*3600;Sl=a*b*pi/T;symse2tf=sqrt(1-e2*cos(t)人2);ft=subs(f,e2,el*el);S=int(ft,0,pi/2);L=4*a*double(S)V=L/TVmax=2*S1/(R+h)Ex22.探月卫星的周期计算程序R=6378;h=[200,600,600,600,600];%输入重要数据%计算长半轴和半焦距%计算离心率和短半轴%定义两个符号变量%定义符号表达式%替换离心率数据%计算积分%符号数据转换为数值%计算平均速度H=[51000,51000,71000,128000,370000];a=(h+H+2*R)/2;c=(H-h)/2;E2=(c./a).人2;b=sqrt(a.*a-c.*c);Vmax=linspace(10.3,10.9,5);bar(Vmax,'c')S=a.*b.*pi;Times=2*S./Vmax./(R+h)/3600Ex23.探月卫星轨道模拟程序functionsatel1()R=6378;h=[200,600,600,600,600];H=[51000,51000,71000,128000,370000];a=(h+H+2*R)/2;c=(H-h)/2;b=sqrt(a.*a-c.*c);t=linspace(-pi,pi,500);x=a'*cos(t);y=b'*sin(t);s=c-c(1);s(2)=-s(2);S=s'*ones(size(t));x=S+x;X=x';Y=y';figure(1),plot(X,Y,'k',-c(1),0,'ob')X=reshape(X,2500,1);Y=reshape(Y,2500,1);figure(2),comet(X,Y)ex24.小行星轨道模拟程序x=[4.5596;5.0816;5.5546;5.9636;6.2756];y=[0.8145;1.3685;1.9895;2.6925;3.5265];A=[x.A2,2*x.*y,y.A2,2*x,2*y];b=-[1;1;1;1;1];ak=A\b;C=[ak(1),ak(2);ak(2),ak(3)];X=-C\[ak(4);ak(5)];x0=X(1);y0=X(2);X=[X;1];D=[ak(1),ak(2),ak(4);ak(2),ak(3),ak(5);ak(4),ak(5),1];F=X'*D*X;[Ud]=eig(C);a=sqrt(-F/d(1,1));b=sqrt(-F/d(2,2));t=2*pi*(0:5000)/5000;u=a*cos(t);v=b*sin(t);V=U*[u;v];x1=V(1,:)+x0;y1=V(2,:)+y0;plot(x1,y1,x,y,'*',x0,y0,'rO'),holdonx2=[x1,x1,x1];y2=[y1,y1,y1];comet(x2,y2)disp([x0,y0])disp([a,b])ex25.汽车租赁业务的矩阵特征向量应用程序X=[120;150];p=0.1;q=0.12;A=[1-p,q;p1-q];Cars=X;fork=1:6X=A*X;Cars=[Cars,X];endCarsfigure(1),bar(Cars(1,:))figure(2),bar(Cars(2,:))alpha=[q;p];R=alpha/sum(alpha);X0=R*270;Taxs=X0;fork=1:6X0=A*X0;Taxs=[Taxs,X0];endTaxsfigure(3),bar(Taxs(1,:))figure(4),bar(Taxs(2,:))ex26.汽车紧急刹车滑行距离的数据拟合R1=1.609;R2=0.3048;x=[2025303540455055606570];

y=[202841537293118149y=[202841537293118149P2=polyfit(x,y,2);

y2=polyval(P2,x);Rest2=sum((y-y2).人2)figure(1),plot(x,y,'*',x,y2)P3=polyfit(x,y,3);y3=polyval(P3,x);Rest3=sum((y-y3).人2)figure(2),plot(x,y,'ok',x,y3,'k')V=[20,30,40,60,80,100,120]./Rl;S=polyval(P3,V);[x;y][V*R1;S*R2]ex27.两船相遇的随机模拟实验程序1266];%二次多项式拟合%计算二次拟合多项式在数据点处的值%计算二次拟合多项式的残差平方和%三次多项式拟合%计算三次拟合多项式在数据点处的值%计算三次拟合多项式的残差平方和%将公里数转换为英里数%使用三次多项式计算刹车距离%显示公里数和刹车距离function[P,Fn]=shipmeet(N)ifnargin==0,N=10000;endS1=0.5*22*22;S2=0.5*23*23;P=1-(S1+S2)/(24*24);data=24*rand(N,2);X=data(:,1);Y=data(:,2);II=find(X-1<=Y&Y<=X+2);Fn=length(II)/N;plot(X(II),Y(II),'.')ex28.平面多边形填充和定积分蒙特卡罗计算实验x1=0:.01:1;y1=sqrt(x1);x2=1:-.01:0;y2=x2.A2;fill([x1,x2],[y1,y2],'r')N=10000;fork=1:6data=rand(N,2);x=data(:,1);y=data(:,2);II=find(y<=sqrt(x)&y>=x.A2);M=length(II);S(k)=M/N;endSsymsuf=sqrt(u)-uA2;Si=int(f,u,0,1)Ex29.维维安尼绘图及体积计算程序r=linspace(0,1,20);alfa=linspace(0,2*pi,30);X=1+r'*cos(alfa);Y=r'*sin(alfa);Z=sqrt(4.01-X.A2-Y.A2);meshz(X,Y,Z),axisoffcolormap([0,0,1])view(133,16)symsxy;f=sqrt(4-xA2-yA2);y1=-sqrt(2*x-xA2);y2=sqrt(2*x-xA2);S1=int(f,y,y1,y2);S2=int(S1,x,0,2)Ex30.维维安尼体积计算蒙特卡罗实验function[V,error]=MonteC(L)ifnargin==0,L=7;endfork=1:LS(k)=viviani;endV=mean(S);error=V-(8/3*pi-32/9);hist(S,7)functionV=viviani(n)ifnargin==0,n=10000;endP=rand(n,3);X=2*P(:,1);Y=2*P(:,2)-1;Z=2*P(:,3);II=find((X-1).A2+Y.A2<=1&Z<=sqrt(4-X.A2-Y.A2));V=8*length(II)/n;Ex31.维维安尼体动态积切割演示程序functionvivia0[X,Y,Z]=sphere(36);colormap([000])mesh(X,Y,Z),axisoffview(150,24)holdont=linspace(0,2*pi,36);xt=.5*cos(t)+.5;yt=.5*sin(t);x=[xt;xt];y=[yt;yt];z=[-1;1.1]*ones(1,36);II=31:36;fork=30:-1:1II=[k,II];mesh(x(:,II),y(:,II),z(:,II)),pause(.5)endzt=sqrt(l+eps-xt42-yt.A2);plot3([xt,xt],[yt,yt],[zt,-zt],'k')ex32.贷款与还贷问题模拟程序fun=inline('100*1.1A(x+1)-300*(1.1Ax-1)');x0=fzero(fun,5);x0+1p=1.1;S=100;k=1;S=S*p;pay=S;whileS>0k=k+1;S=S*p;S=S-30;pay=[pay,S];endpayex33.梯子问题的数学模型求解程序fun=inline('3./sin(alpha)+2./cos(alpha)');figure(1)fplot(fun,[pi/6,2*pi/5]),holdon[x0,y0]=fminbnd(fun,0,pi/2)plot(x0,y0,'o')fun1=inline('(1+2./x).*sqrt(x.A2+9)');figure(2)fplot(fun1,[2,4]),holdon[x1,y1]=fminbnd(fun1,2,4)plot(x1,y1,'o')ex34.单位根实验functionR=myroots(n)ifnargin==0,n=3;endt=linspace(0,2*pi,40);x=cos(t);y=sin(t);P(1)=1;P(n+1)=1;R=roots(P);X=real(R);Y=imag(R);plot(x,y,X,Y,'ro')ex35•常微分方程向量场模拟程序[x,y]=meshgrid(0:.25:6,0:.05:2);k=y.*(1-y);d=sqrt(1+k.A2);px=1./d;py=k./d;quiver(x,y,px,py),holdonu=dsolve('Du=u*(1-u)','u(0)=.2');v=dsolve('Dv=v*(1-v)','v(0)=1.8');ezplot(u,[0,6])ezplot(v,[0,6])axis([0,6.3,-.04,2.01])ex36.追击问题的模拟实验程序function[T,Lp,Lq]=chase(dt)ifnargin==0,dt=1;endV1=1;V2=2*V1;P=[100,0];Q=[0,0];e=[-1,0];d=100;k=1;whiled>0.5P(k+1,:)=P(k,:)+V2*dt*e;Q(k+1,2)=Q(k,2)+V1*dt;e=Q(k+1,:)-P(k+1,:);d=norm(e);e=e/d;k=k+1;endu=Q(:,1);v=Q(:,2);x=P(:,1);y=P(:,2);plot(x,y,'ro',u,v,'b.')T=k*dt;Lp=k*V2*dt;Lq=k*V1*dt;Ex37.渡船航线模拟实验程序function[distance,times]=route(v)ifnargin==0,v=0.5;endB=[0,100];V=1;dt=1;x=v;y=V;distance=sqrt(x^2+(100-y)人2);P=[x,y];times=1;whiledistance>0.5x=x+dt*(v-V