当前位置：文档库 › Matlab在物理上的应用举例

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+z^2); %

Phi=2*pi.*(L-z)./Lambda; %

SumCos=sum(cos(Phi)); %

SumSin=sum(sin(Phi)); %

B(j)=(SumCos^2+SumSin^2)/NPoints^2; %

end

clf,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);

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); end

end

clf;

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 all

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(-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);

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

end

Bax=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 对电荷#%g\n',ic);

rc=input('电荷位置[x y]（米）：');

x(ic)=rc(1); %

y(ic)=rc(2); %

q(ic)=input('输入电荷量（库仑）'); end

E0=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(xij^2+yij^2); %

Fx=Fx+C0*q(ic)*xij/Rij^3;

Fy=Fy+C0*q(ic)*yij/Rij^3;

end

fprintf('其它电荷作用在电荷#%g上的合力为：\n',ic); %

fprintf('x－分量：%gN\n',Fx);

fprintf('y－分量：%gN\n',Fy)

end

７．李萨如图形

% lisaru.m

syms t a1 a2 w1 w2

x=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个点,求出时间增量dt

xplot=eval(x);yplot=eval(y); %计算Ns个点的位置x(t),y(t)

figure(gcf);

subplot(1,2,1),

for i=1:750

plot(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=A\B;

a1=X(1), a2=X(2), N1=X(3),N2=X(4)

１１。循环

%%xunhuan

R=8.31;

gama=1.4; %(注:由于在MA TLAB中无拉丁文,所以用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），%气体常数R＝8.314；给定起始总功WTotal=0;点序号iPoint=1;画等温线用

%的点数NCurve=100;P－V图第一点坐标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('对过程＃％g\n',iPoint-1);

PathType=menu(sprintf('过程％g：选择下一路径',iPoint-1),...

'等压','等容','等温','绝热','退出'); %图形界面菜单生成语句

switch PathType

case 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(Q>0)

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

%上两句用于美观图象，使曲线不要紧靠边框

%标注语句略

end

WTotal

QTotal

y=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:2

alfa=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('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.弹簧振子的物理模型分析

本题目用MA TLAB处理简单弹簧振子和耦合弹簧振子的物理模型，并对实验数据进行了拟合。

(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）耦合弹簧振子

用劲度系数为k1的边弹簧和劲度系数为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

L=U’

②在原始方程M*x’’+K*x=0两端各左乘以U’及在中间的对角矩阵U’*U，得

U’*M* U’*U*x’’+ U’*K*U *x=0

作变量置换

z=[z1,z2]t= U’*[x1,x2] t= U’*x,

得

z’’+L*z=0

这是一个对角矩阵方程，既可把它分为两个方程

z1’’+L*z1=0

z2’’+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,xd02]t，x0=[x01,x02]t，则最后得到结果:.

x(t)=[u1]*([u1]t*M*x0*cos(ω1*t)+1/w1*[u1]t*M*xd0*sin(ω1*t))

[u2] *([u2] t*M*x0*cos(ω1*t)+1/w2*[u2]t*M*xd0*sin(ω2*t))

几种不同情况的MA TLAB编程

①同方向的振动,在t=0时刻,两振子从平衡位置向同一方向移动相同的位移,对此质量相同的两振子来说,两振子以同样的频率,同方向等振幅振动.此时两振子的距离保持常数没有能量的传递,两振子的振动频率仅与两边弹簧有关.即两振子同向振动，此时耦合弹簧不起作用.

已知条件：两振子的质量: m1=2;m2=2;

边弹簧的劲度系数：K1=8；

耦合弹簧的劲度系数：K2=4;

两振子的初始位置（在x轴上）：m1的初始位置：x01=1；

m2的初始位置：x02=1

两振子的初始速度（以x轴正向取正）：m1的初始速度：x0d1=6；

m2的初始速度：x0d2=6；

MA TLAB编程

clear all

m1=.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:2

alfa=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

for r=1:2 %分别对x1,x1绘图

subplot(2,1,r)

plot(t,x(r,:)),grid;

xlabel('时间');

ylabel(['位移',num2str(r)]); %绘图

end

运行结果：

w =

2.0000 2.8284

w =

2.0000 2.8284

w =

2.0000 2.8284

w = 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 all

m1=.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); %时间分割和输出变量初始化f or s=1:2 %分别处理两特征值

for j=1:2

alfa=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

for r=1:2 %分别对x1,x1绘图subplot(2,1,r)

plot(t,x(r,:)),grid;

xlabel('时间');

ylabel(['位移',num2str(r)]); %绘图

end

运行结果

w =

2.0000 2.8284

w =

2.0000 2.8284

w =

2.0000 2.8284

w =

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 all

m1=.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:2

alfa=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

for 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:2

alfa=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

for 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*s

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,

不同质量的振子振动的周期分别为单位(毫秒)

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]，振子周期的平方为[A T,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=A\B, 有

A=[detax.ones(N,1)]; B=detay;q=A\B

编程

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=A\B;r=1/q(1); %线性拟合

plot(detax,detay,'o'),hold on %绘出原始数据图

xi=0:50;yi=q(1)*xi+q(2);

A1=detax;q0=A1\B; %通过线性拟合

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.m

function 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.fig

gui_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(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});

end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});

else

gui_mainfcn(gui_State, varargin{:});

end

function g_gui9_OpeningFcn(hObject, eventdata, handles, varargin) handles.output = hObject;

guidata(hObject, handles);

function varargout = g_gui9_OutputFcn(hObject, eventdata, handles) varargout{1} = handles.output;

function popupmenu1_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor')); end

function 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', '请输入新体积');

end

function edit1_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor')); end

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

end

data = 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')); end

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

end

data = 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')); end

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

end

data = getappdata(gcbf, 'metricdata');

data.edit3 = edit3;

setappdata(gcbf, 'metricdata', data);

function edit4_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor')); end

function edit4 = edit4_Callback(hObject, eventdata, handles)

edit4 = str2double(get(hObject, 'String'));

if isnan(edit4)|edit4==0

set(hObject, 'String', '');

errordlg('Input must be a number but 0','Error');

end

data = getappdata(gcbf, 'metricdata');

data.edit4 = edit4;

setappdata(gcbf, 'metricdata', data);

function edit5_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');

else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));

end

function edit5 = edit5_Callback(hObject, eventdata, handles)

edit5 = str2double(get(hObject, 'String'));

if isnan(edit5)|edit5==0

set(hObject, 'String', '');

errordlg('Input must be a number but 0','Error');

end

data = getappdata(gcbf, 'metricdata');

data.edit5 = edit5;

setappdata(gcbf, 'metricdata', data);

function pushbutton1_Callback(hObject, eventdata, handles)

data = getappdata(gcbf, 'metricdata');

popup_sel_index = get(handles.popupmenu1, 'Value');

else

switch popup_sel_index

case 1 % Choice

warndlg('请选择合适的路径','Warning');

case 2 % Deng Ya

if data.edit2==data.edit5

warndlg('请输入新体积','Warning');

else

edit2=min(data.edit2,data.edit5);

edit5=max(data.edit2,data.edit5);

plot([edit2 edit5],[data.edit1 data.edit1]);

grid on;

hold on;

edit3=data.edit3./data.edit2.*data.edit5;

data.edit7=data.edit7+data.edit1.*(data.edit2-data.edit5);

edit6=data.edit1.*data.edit2./data.edit3.*(edit3-data.edit3).*3.5;

data.edit6=data.edit6+edit6;

data.edit8=data.edit8+(edit6+abs(edit6))./2;

data.edit9=data.edit9-(edit6-abs(edit6))./2;

data.edit3=edit3;

data.edit2=data.edit5;

data.edit4=data.edit1;

setappdata(gcbf, 'metricdata', data);

set(handles.edit1, 'String', data.edit1);

set(handles.edit2, 'String', data.edit2);

set(handles.edit3, 'String', data.edit3);

set(handles.edit4, 'String', data.edit4);

set(handles.edit5, 'String', data.edit5);

set(handles.edit6, 'String', data.edit6);

set(handles.edit7, 'String', data.edit7);

set(handles.edit8, 'String', data.edit8);

set(handles.edit9, 'String', data.edit9);

end

case 3 % Deng Rong

if data.edit1==data.edit4

warndlg('请输入新压强','Warning');

else

edit1=min(data.edit1,data.edit4);

edit4=max(data.edit1,data.edit4);

plot([data.edit2 data.edit2],[edit1 edit4]);

grid on;

hold on;

edit3=data.edit3./data.edit1.*data.edit4;

edit6=data.edit1.*data.edit2./data.edit3.*(edit3-data.edit3).*2.5;

data.edit6=data.edit6+edit6;

data.edit8=data.edit8+(edit6+abs(edit6))./2;

data.edit9=data.edit9-(edit6-abs(edit6))./2;

data.edit3=edit3;

data.edit1=data.edit4;

data.edit5=data.edit2;

setappdata(gcbf, 'metricdata', data);

set(handles.edit1, 'String', data.edit1);

set(handles.edit2, 'String', data.edit2);

set(handles.edit3, 'String', data.edit3);

set(handles.edit4, 'String', data.edit4);

set(handles.edit5, 'String', data.edit5);

set(handles.edit6, 'String', data.edit6);

set(handles.edit7, 'String', data.edit7);

set(handles.edit8, 'String', data.edit8);

set(handles.edit9, 'String', data.edit9);

end

case 4 % Deng Wen

if data.edit2==data.edit5

warndlg('请输入新体积','Warning');

else

edit2=min(data.edit2,data.edit5);

edit5=max(data.edit2,data.edit5);

jiange=edit5./1000;

x=[edit2:jiange:edit5];

y=data.edit1.*data.edit2./x;

plot(x,y);

grid on;

hold on;

edit6=data.edit1.*data.edit2.*log(data.edit5./data.edit2);

data.edit6=data.edit6+edit6;

data.edit7=data.edit7-data.edit1.*data.edit2.*log(data.edit5./data.edit2);

data.edit8=data.edit8+(edit6+abs(edit6))./2;

data.edit9=data.edit9-(edit6-abs(edit6))./2;

data.edit4=data.edit1.*data.edit2./data.edit5;

data.edit1=data.edit4;

data.edit2=data.edit5;

setappdata(gcbf, 'metricdata', data);

set(handles.edit1, 'String', data.edit1);

set(handles.edit2, 'String', data.edit2);

set(handles.edit3, 'String', data.edit3);

set(handles.edit4, 'String', data.edit4);

set(handles.edit5, 'String', data.edit5);

set(handles.edit6, 'String', data.edit6);

set(handles.edit7, 'String', data.edit7);

set(handles.edit8, 'String', data.edit8);

set(handles.edit9, 'String', data.edit9);

end

case 5 % Jue Re

if data.edit2==data.edit5

warndlg('请输入新体积','Warning');

else

edit2=min(data.edit2,data.edit5);

edit5=max(data.edit2,data.edit5);

jiange=edit5./1000;

x=[edit2:jiange:edit5];

y=data.edit1./(x.^1.4).*(data.edit2.^1.4);

plot(x,y);

grid on;

hold on;

data.edit4=data.edit1.*data.edit2.^1.4./data.edit5.^1.4;