Matlab在物理上的应用举例

1. 单列波%单列波t=0:0.001:10;A=input(振幅A=);w=input(频率w=);a=input(a=);y=A.*sin(w.*t+a);plot(t,y);pause(1),sound(y);ylabel(y),xlabel(t)2. %光的单缝衍射现象Lambda=500e-9; % a=input(a=); % 可取0.2e-3,1e-3,2e-3三种情况z=1 %ymax=3*Lambda*z/a; %Ny=51; %ys=linspace(-ymax,ymax,Ny); %NPoints=51; %yPoint=linspace(-a/2,a/2,NPoints); % for j=1:Ny % L=sqrt(ys(j)-yPoint).2+z2); % Phi=2*pi.*(L-z)./Lambda; % SumCos=sum(cos(Phi); % SumSin=sum(sin(Phi); % B(j)=(SumCos2+SumSin2)/NPoints2; %endclf,plot(ys,B,*,ys,B);grid; %3. %用毕奥沙伐尔定律计算电流环产生的磁场mu0=4*pi*1e-7;I0=5.0;Rh=1;C0=mu0/(4*pi)*I0;NGx=21;NGy=21;x=linspace(-Rh,Rh,NGx);y=linspace(-3,3,20);y=x;Nh=20;theta0=linspace(0,2*pi,Nh+1);theta1=theta0(1:Nh);y1=Rh*cos(theta1);z1=Rh*sin(theta1);theta2=theta0(2:Nh+1);y2=Rh*cos(theta2);z2=Rh*sin(theta2);dlx=0;dly=y2-y1;dlz=z2-z1;xc=0;yc=(y2+y1)/2;zc=(z2+z1)/2;for i=1:NGy for j=1:NGx rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; r3=sqrt(rx.2+ry.2+rz.2).3; dlXr_x=dly.*rz-dlz.*ry; dlXr_y=dly.*rx-dlx.*rz; Bx(i,j)=sum(C0*dlXr_x./r3); By(i,j)=sum(C0*dlXr_y./r3);endendclf;quiver(x,y,Bx,By);4. %多普勒效应x0=500;v=50;y0=20; c=330;w=1000; t=0:0.001:30; r=sqrt(x0-v*t).2+y0.2); t1=t-r/c; u=sin(w*t)+sin(1.1*w*t); u1=sin(w*t1)+sin(1.1*w*t1); sound(u);pause(5);sound(u1);5.亥姆霍兹线圈 clear allmu0=4*pi*1e-7;I0=5.0;Rh=1;C0=mu0/(4*pi)*I0;NGx=21;NGy=21;x=linspace(-Rh,Rh,NGx);y=linspace(-Rh,Rh,NGy);Nh=20;theta0=linspace(0,2*pi,Nh+1);theta1=theta0(1:Nh);y1=Rh*cos(theta1);z1=Rh*sin(theta1);theta2=theta0(2:Nh+1);y2=Rh*cos(theta2);z2=Rh*sin(theta2);dlx=0;dly=y2-y1;dlz=z2-z1;xc=0;yc=(y2+y1)/2;zc=(z2+z1)/2;for i=1:NGy for j=1:NGx rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; r3=sqrt(rx.2+ry.2+rz.2).3; dlXr_x=dly.*rz-dlz.*ry; dlXr_y=dly.*rx-dlx.*rz; Bx(i,j)=sum(C0*dlXr_x./r3); By(i,j)=sum(C0*dlXr_y./r3);endendBax=Bx(:,11:21)+Bx(:,1:11);Bay=By(:,11:21)+By(:,1:11);subplot(1,2,1)mesh(x(11:21),y,Bax);xlabel(x);ylabel(y);subplot(1,2,2),plot(y,Bax),grid,xlabel(y);ylabel(Bx);库仑引力clear all;N=input(电荷数目N：);for ic=1:N % fprintf(-n 对电荷#%gn,ic); rc=input(电荷位置xy（米）：); x(ic)=rc(1); % y(ic)=rc(2); % q(ic)=input(输入电荷量（库仑）);endE0=8.85e-12; %C0=1/(4*pi*E0); %for ic=1:N Fx=0.0; Fy=0:0; % for jc=1:N % if(ic=jc) % xij=x(ic)-x(jc); yij=y(ic)-y(jc); Rij=sqrt(xij2+yij2); % Fx=Fx+C0*q(ic)*xij/Rij3; Fy=Fy+C0*q(ic)*yij/Rij3; end end fprintf(其它电荷作用在电荷#%g上的合力为：n,ic); % fprintf(x分量：%gNn,Fx); fprintf(y分量：%gNn,Fy)end李萨如图形 % lisaru.msyms t a1 a2 w1 w2x=cos(w1.*t+a1);y=sin(w2.*t+a2);a1=input(a1=);a2=input(a2=);w1=input(w1=);w2=input(w2=);tf=10;Ns=1000;t=linspace(0,tf,Ns); dt=tf/(Ns-1); %分Ns个点,求出时间增量dtxplot=eval(x);yplot=eval(y); %计算Ns个点的位置x(t),y(t)figure(gcf);subplot(1,2,1),for i=1:750plot(yplot(1:i),xplot(1:i); %画点的轨迹图axis(equal); grid ; %使两轴比例相同pause(0.01)end耦合振子 m1=2;m2=2;K1=16;K2=4; %?x0=1;0;xd0=6;6;tf=10; %?M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ; %?u(:,s)u,L=eig(K,M) ; %?t=linspace(0,tf,101);x=zeros(2,101); %?for s=1:2 %?alfa=sqrt(u(:,s)*M*u(:,s) ; %?u(:,s)= u(:,s)/alfa; %? w(j)=sqrt(L(j,j); %?xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt; %?end for r=1:2 %?x1,x1?subplot(2,1,r)plot(t.x(r,:),grid;xlabel(xxx);ylabel(yyy,num2str(r);end 拍频%t=0:0.001:10; a1=input(?1=);w1=input(?1=); a2=input(?2=);w2=input(?2=);y1=a1*sin(w1*t); y2=a2*sin(w2*t);y=y1+y2; subplot(3,1,1),plot(t,y1),ylabel(y1) subplot(3,1,2),plot(t,y2),ylabel(y2)subplot(3,1,3),plot(t,y),ylabel(y),xlabel(t)pause,sound(y1);pause(5),sound(y2);pause(5),sound(y),pause subplot(1,1,1) 。物块下滑m1=input(m1=);m2=input(m2=);theta=input(theta(度)=);theta=theta*pi/180; g=9.81;A=m1*cos(theta),-m1,-sin(theta),0;. m1*sin(theta),0,cos(theta),0;. 0,m2, -sin(theta),0;. 0,0,-cos(theta),1;B=0,m1*g,0,m2*g; X=AB;a1=X(1), a2=X(2), N1=X(3),N2=X(4) 。循环%xunhuanR=8.31;gama=1.4; %(注:由于在MATLAB中无拉丁文,所以用gama代替?) nMoles=0.5;P(1)=1e5;V(1)=0.012;WTotal=0;QTotal=0;iPoint=1;NCurve=100;PPlot=P(1);VPlot=V(1);%变量和图形初始化，输入气体的摩尔数nMoles，初始压力P（1），初始容积V（1），%气体常数R8.314；给定起始总功WTotal=0;点序号iPoint=1;画等温线用%的点数NCurve=100;PV图第一点坐标PPlot=P(1);VPlot=V(1)T(1)=P(1)*V(1)/(nMoles*R); %算出初始温度%为了进入循环，先要设两个不相等的PathType和QuitType值QuitType=5;PathType=0;while(PathType=QuitType) %在菜单上选择退出之前不断循环， %选择路径类型或退出 iPoint=iPoint+1; %下一点 fprintf(对过程gn,iPoint-1); PathType=menu(sprintf(过程g：选择下一路径,iPoint-1),. 等压,等容,等温,绝热,退出); %图形界面菜单生成语句 switch PathTypecase 1 %等压路径 V(iPoint)=input(输入新容积：); P(iPoint)=P(iPoint-1); %压力不变 T(iPoint)=P(iPoint)*V(iPoint)/(nMoles*R); %按新容积算出温度 W=P(iPoint)*(V(iPoint)-V(iPoint-1); %计算等压过程所做的功 Q=(gama*nMoles*R/(gama-1)*(T(iPoint)-T(iPoint-1);VPlot=VPlot,V(iPoint); %加上新的容积和压力点，用以绘图 PPlot=PPlot,P(iPoint);case 2 %等容路径 P(iPoint)=input(输入新压力：); V(iPoint)=V(iPoint-1); %容积不变 T(iPoint)=P(iPoint)*V(iPoint)/(nMoles*R); %按新压力算出温度 W=0; %等容路径上所做的功为零 Q=(nMoles*R/(gama-1)*(T(iPoint)-T(iPoint-1); VPlot=VPlot,V(iPoint); %加上绘图用的新容积和压力点 PPlot=PPlot,P(iPoint);case 3 %等温路径 V(iPoint)=input(输入新容积：); T(iPoint)=T(iPoint-1); %温度不变 P(iPoint)=nMoles*R*T(iPoint)/V(iPoint); %按新容积求新压力 W=nMoles*R*T(iPoint)*log(V(iPoint)/V(iPoint-1); %求所做的功 Q=W; %用元素群运算求等温路径上的P和V，加进绘图数据中 VNew=linspace(V(iPoint-1),V(iPoint),NCurve); PNew=nMoles*R*T(iPoint)./VNew; VPlot=VPlot,VNew; %将新的V，P点加入绘图数据中 PPlot=PPlot,PNew;case 4 %绝热路径 V(iPoint)=input(输入新容积:); P(iPoint)=P(iPoint-1)*V(iPoint-1)gama/V(iPoint)gama; %按新容积求新压力 T(iPoint)=T(iPoint-1)*V(iPoint-1)(gama-1)/V(iPoint)(gama-1); %按新容积求新温度 W=-(nMoles*R/(gama-1)*(T(iPoint)-T(iPoint-1); %求所做的功 Q=0; VNew=linspace(V(iPoint-1),V(iPoint),NCurve); PNew=P(iPoint-1)*V(iPoint-1)gama./VNew.gama; VPlot=VPlot,VNew; %将新的V，P点加入绘图数据中 PPlot=PPlot,PNew; otherwise end %画出到目前为止的PV图 if(PathType=QuitType) WTotal=WTotal+W; %将新做的功加进总功 if(Q0) QTotal=QTotal+Q; end figure(gcf);plot(V,P,o,VPlot,PPlot,-)%图形窗移前，绘图 size=axis; axis(size(1)*0.9,size(2)*1.1,size(3)*0.9,size(4)*1.1); %上两句用于美观图象，使曲线不要紧靠边框 %标注语句略 endendWTotalQTotaly=WTotal/QTotal 。振子和拍m1=2;m2=2;K1=8;K2=4; %输入各原始参数x0=1;0;xd0=6;0;tf=10; %初使条件M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ; %构成参数矩阵u(:,s)*u,L=eig(K,M) ; %求广义特征向量和特征值t=linspace(0,tf,101);x=zeros(2,101); %时间分割和输出变量初始化for s=1:2 %分别处理两特征值 for j=1:2alfa=sqrt(u(:,s)*M*u(:,s) ; %解耦后的向量u(:,s)= u(:,s)/alfa; %把特征向量归一化 m=sqrt(L); %分别求对应于两特征值的分量w(j)=m(j,j)xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt; %把两个分量累加endendfor r=1:2 %分别对x1,x1绘图subplot(2,1,r)plot(t,x(r,:),grid;xlabel(xxx);ylabel(yyy,num2str(r);end 。阻尼振动B=0.7;w=5;t=0:0.001:10;x=dsolve(D2x)+2.*B.*(Dx)+w.2.*x);plot(t,x)7.弹簧振子的物理模型分析l 本题目用MATLAB处理简单弹簧振子和耦合弹簧振子的物理模型，并对实验数据进行了拟合。n(1) 简单弹簧振子用劲度系数分别为k1和k2的弹簧，把质量为m的振子固定在气垫导轨，求振子的速度v,加速度a和位移x。根据牛顿第二定律m*(D2x)=-(k1+k2)*x解微分方程，得运动方程x=A*cos(w*t+)对上式微分得速度和加速度v=-A*w*sin(w*t+)a=-A*w*w*cos(w*t+)其中 w=sqrt(k1+k2)/m); 编程 syms x m k1 k2;x=dsolve (m*(D2x)+(k1+k2)*x=0,x(0)=0,Dx(0)=10,t);v=diff(x)a=diff(v)（2）耦合弹簧振子l 用劲度系数为k 1的边弹簧和劲度系数为k2的耦合弹簧连接质量为m1和m2的两振子，求此物理模态的振动系统方程。设x1和x2分别是两振子离开平衡位置的位移，耦合振子的振动方程是m1*(d(dx1/dt)/dt)=-k1*x1-k2*(x1-x2);m2*(d(dx2/dt)/dt)=-k1*x2-k2*(x2-x1);写成矩阵形式M*(D2x)+K*x=0 用eig函数求出矩阵K-M的特征值L和特征向量U, U和L满足U*M*U=I 1 00 2L=U*K*U= 在原始方程M*x+K*x=0两端各左乘以U及在中间的对角矩阵U*U，得U*M* U*U*x+ U*K*U *x=0作变量置换z=z1,z2t= U*x1,x2 t= U*x,得z+L*z=0这是一个对角矩阵方程，既可把它分为两个方程z1+L*z1=0z2+L*z2=0这意味着两种振动模态可以解耦.1*1=1*1，2*2=2*2，其中1是第1个模态的固有频率,2是第2个模态的固有频率. 由上述的解耦模态中，给出初始条件x0,xd0，化为z0,zd0求出变量z1,z2再把z转换为x. 设速度和位置的初始条件分别为xd0=xd01,xd02t， x0=x01,x02t，则最后得到结果:.x(t)=u1*(u1t*M*x0*cos(1*t)+1/w1*u1t*M*xd0*sin(1*t)u2 *(u2 t*M*x0*cos(1*t)+1/w2*u2t*M*xd0*sin(2*t) 几种不同情况的MATLAB编程 同方向的振动,在t=0时刻,两振子从平衡位置向同一方向移动相同的位移,对此质量相同的两振子来说,两振子以同样的频率,同方向等振幅振动.此时两振子的距离保持常数没有能量的传递,两振子的振动频率仅与两边弹簧有关.即两振子同向振动，此时耦合弹簧不起作用.已知条件： 两振子的质量: m1=2;m2=2;边弹簧的劲度系数：K1=8； 耦合弹簧的劲度系数：K2=4; 两振子的初始位置（在x轴上）：m1的初始位置：x01=1；m2的初始位置：x02=1两振子的初始速度（以x轴正向取正）：m1的初始速度：x0d1=6；m2的初始速度：x0d2=6；MATLAB编程 clear allm1=.2;m2=.2;K1=8;K2=8; %输入各原始参数x0=1;0;xd0=6;6;tf=10; %初始条件M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ; %构成参数矩阵u,L=eig(K,M) ; %求广义特征向量和特征值t=linspace(0,tf,101);x=zeros(2,101); %时间分割和输出变量初始化for s=1:2 %分别处理两特征值 for j=1:2alfa=sqrt(u(:,s)*M*u(:,s) ; %解耦后的向量u(:,s)= u(:,s)/alfa; %把特征向量归一化 m=sqrt(L); %分别求对应于两特征值的分量w(j)=m(j,j)xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt; %把两个分量累加end end for r=1:2 %分别对x1,x1绘图subplot(2,1,r)plot(t,x(r,:),grid;xlabel(时间);ylabel(位移,num2str(r); %绘图end 运行结果：w = 2.0000 2.8284w = 2.0000 2.8284w = 2.0000 2.8284w = 2.0000 2.8284 振子反方向振动时,在t=0时,使质量相同的两振子位于各自的平衡位置两侧,且到平衡位置的距离相等,此后两振子以相同的频率,方向等幅振动.两振子不发生能量传递.但此时的频率比同方向振动的频率大,周期要短,因为此时耦合弹簧起作用.已知条件 ： 两振子的质量: m1=2;m2=2; 边弹簧的劲度系数：K1=8；耦合弹簧的劲度系数：K2=4; 两振子的初始位置（在x轴上建坐标系）：m1的初始位置：x01=1；m2的初始位置：x02=1；两振子的初始速度（以x轴正向取正值）：m1的初始速度：x0d1=6；m2的初始速度：x0d2=-6；编程clear allm1=.2;m2=.2;K1=4;K2=4; %输入各原始参数x0=1;0;xd0=6;-6;tf=40; %初使条件M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ; %构成参数矩阵u( u,L=eig(K,M) ; %求广义特征向量和特征值t=linspace(0,tf,101);x=zeros(2,101); %时间分割和输出变量初始化for s=1:2 %分别处理两特征值 for j=1:2alfa=sqrt(u(:,s)*M*u(:,s) ; %解耦后的向量u(:,s)= u(:,s)/alfa; %把特征向量归一化 m=sqrt(L); %分别求对应于两特征值的分量w(j)=m(j,j)xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt; %把两个分量累加end endfor r=1:2 %分别对x1,x1绘图subplot(2,1,r)plot(t,x(r,:),grid; xlabel(时间);ylabel(位移,num2str(r); %绘图 end 运行结果w = 2.0000 2.8284w = 2.0000 2.8284w = 2.0000 2.8284w = 2.0000 2.8284 拍振动,在t=0时,两振子到平衡位置的距离不同,振子的质量相同,两振子在耦合弹簧及边弹簧的作用下会发生拍振动.两振子会存在能量传递.当一振子处于平衡位置另一振子离平衡位置距离为A,此两振子在耦合弹簧的作用下会发生拍振动.两振子会存在能量传递.第一个振子振幅变小时,另一个振子振幅增大.已知条件：两振子的质量: m1=2;m2=2;边弹簧的颈度系数：K1=8； 耦合弹簧的颈度系数：K2=4; 两振子的初始位置（在x轴上建坐标系）：m1的初始位置：x01=1；m2的初始位置：x02=1；两振子的初始速度（以x轴正向取正值）：m1的初始速度：x0d1=6；m2的初始速度：x0d2=0；编程clear allm1=.2;m2=.2;K1=8;K2=4; %输入各原始参数x0=1;0;xd0=6;0;tf=40; %初使条件M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ; %构成参数矩阵 u,L=eig(K,M) ; %求广义特征向量和特征值t=linspace(0,tf,101); x=zeros(2,101); %时间分割和输出变量初始化for s=1:2 %分别处理两特征值 for j=1:2alfa=sqrt(u(:,s)*M*u(:,s) ; %解耦后的向量u(:,s)= u(:,s)/alfa; %把特征向量归一化 m=sqrt(L); %分别求对应于两特征值的分量w(j)=m(j,j)xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt; %把两个分量累加end endfor r=1:2 %分别对x1,x1绘图subplot(2,1,r)plot(t,x(r,:),grid; xlabel(时间); ylabel(位移,num2str(r); %绘图 end 在t=0时,使两振子离平衡位置相同,振子的质量不同,两振子在耦合弹簧及边弹簧的作用下，会发生拍振动.两振子之间会存在能量传递已知条件 两振子的质量: m1=2；m2=4；边弹簧的劲度系数：K1=8； 耦合弹簧的劲度系数：K2=4； 两振子的初始位置（在x轴上建坐标系）：m1的初始位置：x01=1；m2的初始位置：x02=1两振子的初始速度（以x轴正向取正值）：m1的初始速度：x0d1=6；m2的初始速度：x0d2=6；编程： m1=.2;m2=.4;K1=8;K2=4; %输入各原始参数x0=1;0;xd0=6;6;tf=40; %初使条件M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ; %构成参数矩阵u,L=eig(K,M) ; %求广义特征向量和特征值t=linspace(0,tf,101);x=zeros(2,101); %时间分割和输出变量初始化for s=1:2 %分别处理两特征值 for j=1:2alfa=sqrt(u(:,s)*M*u(:,s) ; %解耦后的向量u(:,s)= u(:,s)/alfa; %把特征向量归一化 m=sqrt(L); %分别求对应于两特征值的分量w(j)=m(j,j)xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt; %把两个分量累加end endfor r=1:2 %分别对x1,x1绘图subplot(2,1,r)plot(t,x(r,:),grid;xlabel(时间); ylabel(位移,num2str(r); %绘图 end (3)弹簧振子有关物理量测量及实验数据分析实验数据记录 :弹簧弹性系数s=8.64cm, m=21克,g=979.44cm/s*sm0=21.830, m1=187.860, m2=49.790, m3=99.540, m4=99.790,振子的质量分别为克。 a=m1+m0/3,b=m1+m2+m0/3,c=m1+m3+m0/3,d=m1+m2+m3+m0/3,e=m1+m3+m4+m0/3,不同质量的振子振动的周期分别为 单位(毫秒)at=(38446.6+38447.0+38449.2)/30,bt=(43077.1+43078.8+43078.1)/30,ct=(47256.0+47255.8+47256.1)/30,dt=(51078.4+51079.3+51078.3)/30,et=(54670.9+54668.0+54667.6)/30,振子的质量是a,b,c,d,e测得的振子周期为at,bt,ct,dt.et，振子周期的平方为AT,BT,CT,DT,ET,FT又，振子的周期的平方与振子的质量成正比,以下模型利用实验数据拟合二次曲线设直线的方程是y=a(1)*=x+a(2)待定的系数是a(1),a(2).将实验数据分别带入x,y,得方程组，那么设方程组的系数矩阵为detax,detay。detax*q(1)+ones(N,1)*q(2)=detail其中，detax,detay均为列,N个一次方程，方程组中只有两个未知数，是超定方程，用最小二乘法可以直接运算q=AB, 有A=detax.ones(N,1); B=detay;q=AB编程m0=21.830; m1=187.860; m2=49.790; m3=99.540; m4=99.790;a=m1+m0/3; b=m1+m2+m0/3; c=m1+m3+m0/3; d=m1+m2+m3+m0/3; e=m1+m3+m4+m0/3;f=m1+m2+m3+m0/3;at=(38446.6+38447.0+38449.2)/30;bt=(43077.1+43078.8+43078.1)/30;%ct=(47256.0+47255.8+47256.1)/30;dt=(51078.4+51079.3+51078.3)/30;et=(54670.9+54668.0+54667.6)/30;AT=at*at; BT=bt*bt; CT=ct*ct; DT=dt*dt; ET=et*et;detax=AT,BT,CT,DT,ET*10(-6); detay=a,b,c,d,e; %原始数据A=detax,ones(5,1); B=detay;q=AB;r=1/q(1); %线性拟合plot(detax,detay,o),hold on %绘出原始数据图xi=0:50;yi=q(1)*xi+q(2);A1=detax;q0=A1B; %通过线性拟合plot(xi,yi,xi,q0*xi,:) %绘图q2=polyfit(detax,detay,2); yi=polyval(q2,xi) ; %二次拟合plot(xi,yi); xlabel(时间的平方 单位(秒的平方);ylabel(质量 单位(克);hold off 8. 傅科摆文件名：fkb.m%a=60;q=4,0,0,0;c=a*pi/180;t,x=ode45(fkb,0:0.02:220,q, ,c);comet(x(:,1),x(:,3) 文件名：fkbfun.mfunction tt=fkb(t,x,flag,c)a=(2*pi*sin(c)/100;b=9.8/67;tt=x(2); 2*a*x(4)-b*x(1); x(4);-2*a*x(2)-b*x(3);function varargout = g_gui9(varargin)% g_gUI9 M-file for g_gui9.figgui_Singleton = 1;gui_State = struct(gui_Name, mfilename, . gui_Singleton, gui_Singleton, . gui_OpeningFcn, g_gui9_OpeningFcn, . gui_OutputFcn, g_gui9_OutputFcn, . gui_LayoutFcn, , . gui_Callback, );if nargin & isstr(varargin1) gui_State.gui_Callback = str2func(varargin1);endif nargout varargout1:nargout = gui_mainfcn(gui_State, varargin:);else gui_mainfcn(gui_State, varargin:);endfunction g_gui9_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;guidata(hObject, handles);function varargout = g_gui9_OutputFcn(hObject, eventdata, handles)varargout1 = handles.output;function popupmenu1_CreateFcn(hObject, eventdata, handles)if ispc set(hObject,BackgroundColor,white);else set(hObject,BackgroundColor,get(0,defaultUicontrolBackgroundColor);endfunction popupmenu1_Callback(hObject, eventdata, handles)popup_sel_index = get(handles.popupmenu1, Value);switch popup_sel_index case 1 % Choice set(handles.text19, String, 请选择合适的路径); case 2 % Deng Ya set(handles.text19, String, 请输入新体积); case 3 % Deng Rong set(handles.text19, String, 请输入新压强); case 4 % Deng Wen set(handles.text19, String, 请输入新体积); case 5 % Jue Re set(handles.text19, String, 请输入新体积); endfunction edit1_CreateFcn(hObject, eventdata, handles)if ispc set(hObject,BackgroundColor,white);else set(hObject,BackgroundColor,get(0,defaultUicontrolBackgroundColor);endfunction edit1 = edit1_Callback(hObject, eventdata, handles)edit1 = str2double(get(hObject, String);if isnan(edit1)|edit1=0 set(hObject, String, ); errordlg(Input must be a number but 0,Error);enddata = getappdata(gcbf, metricdata);data.edit1 = edit1;setappdata(gcbf, metricdata, data);function edit2_CreateFcn(hObject, eventdata, handles)if ispc set(hObject,BackgroundColor,white);else set(hObject,BackgroundColor,get(0,defaultUicontrolBackgroundColor);endfunction edit2 = edit2_Callback(hObject, eventdata, handles)edit2 = str2double(get(hObject, String);if isnan(edit2)|edit2=0 set(hObject, String, ); errordlg(Input must be a number but 0,Error);enddata = getappdata(gcbf, metricdata);data.edit2 = edit2;setappdata(gcbf, metricdata, data);function edit3_CreateFcn(hObject, eventdata, handles)if ispc set(hObject,BackgroundColor,white);else set(hObject,BackgroundColor,get(0,defaultUicontrolBackgroundColor);endfunction edit3 = edit3_Callback(hObject, eventdata, handles)edit3 = str2double(get(hObject, String);if isnan(edit3)|edit3=0 set(hObject, String, ); errordlg(Input must be a number but 0,Error);enddata = getappdata(gcbf, metricdata);data.edit3 = edit3;setappda