2020iypt，第十三题摩擦振子（二）结果调参代码
代码允许转载，但是放在附录时一定要记得把出处和制作人标清楚！！！！！！！！！！！！！！！！！！！！！！！！！！！
后续的仿真动画代码也会发出。
常见问题：1.代码复制一定要复制到编辑器中，而不是命令行
2.一定保证符号全英文
3.对代码修改后尤其是方程部分修改后，求解错误是正常现象，请一定要首先检查方程是否正确
4.这个代码求解方程是欧拉法求解的方程，如果不了解数值计算的可以先了解一下数值计算！！
使用软件：matlab
代码：
%摩擦振子结果输出动画
%当使用调整变量x：两轴间间距时使用fun1并注释fun2
%当使用调整变量u：两轴间间距时使用fun2并注释fun1
%制作人：茤一分问候 时间2020年/1/14/0：30引用时（这部分请不要删除）
loops=100;
F(loops) = struct('cdata',[],'colormap',[]);
v = VideoWriter('参数u变动.avi');%创建动画
open(v);
for i=0.40:0.0006:0.46
j=round((i-0.4+0.0006)/0.0006);%在改变参数u两轴间间距时使用
%j=round((i-0.2+0.0006)/0.0006);%在改变参数x两轴间间距时使用
figure(j)
fun2(i)%在改变参数u两轴间间距时使用
%fun1(i)%在改变参数x两轴间间距时使用
F(j) = getframe(gcf);%保存帧
writeVideo(v,F(j));%写入动画
close
end
close(v);
function y=fun2(u,y)
%摩擦系数变量调用函数
x=0.20;%两轴的一半间距
l=0.8;%重物的长度
p=1.5*10^3;%重物的密度
d=0.05;%重物的宽度
h=0.01;%重物的高度
g=9.8;%重力加速度
m=p*h*l*d;%重物的质量
%u=0.44;%动摩擦因数
%——————————————————————————
a2=zeros(10000,1);
v2=zeros(10000,1);
x2=zeros(10000,1);
x2(1,1)=0.1;%重心偏移原点的位置
t=0.001;%时间步长
G=m*g;%重力
for i=2:10000
G=m*g;%重力
N1=((x*G-G*x2(i-1,1))/2+x2(i-1,1)*G)/x;%受力分析的方程
N2=(x*G-G*x2(i-1,1))/(2*x); %受力分析的方程
a2(i,1)=(N2*u-N1*u)/m; %受力分析的方程
v2(i,1)=v2(i-1,1)+a2(i,1)*t; %离散速度
x2(i,1)=x2(i-1,1)+v2(i,1)*t; %离散位移
end
hold on
axis([0 10000 -2.5 2.5])
plot(x2,'Color','b') %绘制位移
plot(v2,'Color','r') %绘制速度
plot(a2,'Color','k') %绘制加速度
end
%-------------------------------------------------
function y=fun1(x,y)
%两轴间间距调用变量
%x=0.20;两轴的一半间距
l=0.8;%重物的长度
p=1.5*10^3;%重物的密度
d=0.05;%重物的宽度
h=0.01;%重物的高度
g=9.8;%重力加速度
m=p*h*l*d;%重物的质量
u=0.44;%动摩擦因数
%——————————————————————————
a2=zeros(10000,1);
v2=zeros(10000,1);
x2=zeros(10000,1);
x2(1,1)=0.1;%重心偏移原点的位置
t=0.001;%时间步长
G=m*g;%重力
for i=2:10000
G=m*g;%重力
N1=((x*G-G*x2(i-1,1))/2+x2(i-1,1)*G)/x;%受力分析的方程
N2=(x*G-G*x2(i-1,1))/(2*x); %受力分析的方程
a2(i,1)=(N2*u-N1*u)/m; %受力分析的方程
v2(i,1)=v2(i-1,1)+a2(i,1)*t; %速度
x2(i,1)=x2(i-1,1)+v2(i,1)*t; %位移
end
hold on
axis([0 10000 -2.5 2.5])
plot(x2,'Color','b')
plot(v2,'Color','r')
plot(a2,'Color','k')
end
