大物实验Matlab
黔无驴，有好事者船载以入。
本人没少受各位CSDN大佬、窝工网盘计划各位前辈的恩惠，也想给后来者留点东西。把自己做的几个大物实验的图像处理和数据计算的代码留下吧，万一有人用了呢？
注意：
本文使用的工具为MATLAB R2022a。
—咱不是被ban了吗？怎么搞？
—别问我，问马云（你懂我意思吧，嘿嘿）
451%碰撞打靶试验
2
R=[12.16 14.73 15.92 18.04 18.99 20.83 21.33 23.08 24.46 25.22 27.82 28.87]/100;%水平射程，米
3
y=0.169;%抛出高度，米
4
h0=[5.1 8.5];%打靶高度，厘米
5
H=0.03:0.01:0.14;%高度差，米
6
g=9.8;
7
8
figure(1);%绘制R-h曲线
9
10
plot(100*H,100*R,"-b+",LineWidth=3);%实验散点图绘制
11
hold on;
12
13
h=linspace(2,16,200)/100;%高度差，米
14
r=2.*sqrt(y.*h);%理论射程，米
15
R1=2.*sqrt(y.*H)*100;%各高度理论射程，厘米
16
plot(100*h,100*r,"-r",LineWidth=3);%理论曲线绘制
17
hold on;
18
19
plot(h0,[16 22],"k^",LineWidth=3);
20
grid on;%坐标背景
21
22
xlabel("Δh/cm","FontName","Times New Roman","FontSize",24);
23
ylabel("R(h)/cm","FontName","Times New Roman","FontSize",24);
24
title("图1 高度差h与水平射程R的关系","FontName","宋体","FontSize",24);
25
legend("实际值","理论曲线","打靶点","Location","southeast");
26
hold off;
27
28
figure(2);%绘制碰前碰后速度差曲线
29
30
v=sqrt(2.*g.*H);%碰前速度
31
v1=R*sqrt(g/(2*y));%碰后速度
32
plot(v,v1,"-b*",LineWidth=3);
33
hold on;
34
35
k=polyfit(v,v1,1);%最小二乘法拟合系数
36
v2=linspace(0.7,1.7,100);
37
v3=k(1)*v2+k(2);
38
plot(v2,v3,"-r",LineWidth=3);
39
grid on;
40
41
xlabel("碰前速度差Δv/(m/s)","FontName","宋体","FontSize",24);
42
ylabel("碰后速度差Δv'/(m/s)","FontName","宋体","FontSize",24);
43
title("图2 碰前碰后速度关系","FontName","宋体","FontSize",24);
44
legend("实际值","拟合直线","Location","southeast");
45
hold off;
761%拉伸法测定杨氏弹性模量
2
d_init=[0.627 0.630 0.622 0.624 0.615]/1000;%金属丝直径测量
3
n1=[10.50 10.69 10.88 11.10 11.29 11.45 11.61 11.90 12.08 12.28 12.41]/100;%增重时0~10个0.5kg砝码对应刻度
4
n2=[10.57 10.79 10.91 11.15 11.35 11.41 11.45 11.83 11.90 12.11 12.41]/100;%减重时0~10个0.5kg砝码对应刻度
5
b=8.75/100;
6
l=74/100;
7
L=127/100;
8
delta_b=0.5/1000;
9
delta_L=10/1000;
10
delta_l=5/1000;
11
delta_d=0.004/1000;
12
delta_dif_n=sqrt(0.5)/1000;%以上，米
13
g=9.8066;
14
15
%所需数据计算
16
n=(n1+n2)/2;
17
F=[0 1 2 3 4 5 6 7 8 9 10]*0.5*g;%读数和对应的应力计算
18
d=0;
19
for i=1:5
20
d=d+d_init(i);
21
end
22
d=d/5;%直径平均值，米
23
tp=1.14;
24
25
%逐差
26
sum_del_n=0;
27
sum_del_F=0;
28
del_n=zeros(1,5);
29
del_F=zeros(1,5);
30
for i=1:5
31
del_n(i)=n(i+6)-n(i+1);
32
sum_del_n=sum_del_n+del_n(i);
33
del_F(i)=F(i+6)-F(i+1);
34
sum_del_F=sum_del_F+del_F(i);
35
end%逐差计算delta之和，便于以下操作和输出数据
36
Ey=(8*L*l)/(pi*b*d.^2)*sum_del_F/sum_del_n;%杨氏模量值的计算，Ey便于与相对不确定度区分
37
38
%不确定度准备
39
Sd=std(d_init,1)*tp;%d的A类不确定度
40
ud=delta_d/sqrt(3);%d的B类不确定度
41
Ud=sqrt(Sd.^2+ud.^2);%d的合成不确定度
42
Sn=std(del_n,1)*tp;
43
un=delta_dif_n/sqrt(3);
44
Un=sqrt(Sn.^2++un.^2);%n'-n的合成不确定度
45
Ul=delta_l/sqrt(3);
46
UL=delta_L/sqrt(3);
47
Ub=delta_b/sqrt(3);
48
49
%不确定度计算
50
tmp=L.^(-2)*UL.^2+l.^(-2)*Ul.^2+b.^(-2)*Ub.^2+(2/d).^2*Ud.^2+(5/sum_del_n).^2*Un.^2;
51
E=sqrt(tmp);
52
U=E*Ey;
53
54
%结果输出
55
disp("金属丝直径平均值为（mm）："+d*1000);
56
disp("各次测量标尺读数平均值依次为（cm）：");
57
n*100
58
disp("各次对应应力为（N）：");
59
F
60
disp("逐差各部分读数差依次为（cm）：");
61
del_n*100
62
disp("逐差拉力差总和为"+sum_del_F+"N");
63
disp("逐差读数差总和为"+sum_del_n*100+"cm");
64
disp("测得金属丝杨氏模量为（N/m2）："+Ey);
65
disp("直径的A类不确定度为"+Sd*1000+"mm");
66
disp("直径的B类不确定度为"+ud*1000+"mm");
67
disp("直径的合成不确定度为"+Ud*1000+"mm");
68
disp("n'-n的A类不确定度为"+Sn*100+"cm");
69
disp("n'-n的B类不确定度为"+un*100+"cm");
70
disp("n'-n的合成不确定度为"+Un*100+"cm");
71
disp("金属丝长度l的B类不确定度为"+Ul*100+"cm");
72
disp("镜子与望远镜距离L的B类不确定度为"+UL*100+"cm");
73
disp("镜子到金属丝长度b的B类不确定度为"+Ub*100+"cm");
74
disp("测得杨氏模量的相对不确定度为"+E*100+"%");
75
disp("杨氏模量绝对不确定度为"+U+"N/m2");
76
disp("最终彻底，金属丝杨氏模量表达式为"+Ey+"±"+U+"（N/m2）");
731%液体粘度的测定
2
x1=[25.435 14.810 8.312 20.356 15.000 15.420 20.936 20.595 19.374 28.492 17.029 25.581 32.498 24.885 32.130 21.445 31.764 42.852 32.194 39.482 20.378 29.146 37.975 44.934 41.480];%毫米
3
x2=[26.832 16.238 9.728 21.770 16.432 16.841 22.336 22.005 20.810 29.912 18.442 27.000 33.876 26.324 33.558 22.859 33.159 44.264 33.620 40.899 21.808 30.600 39.462 46.365 42.896];%毫米
4
T=[25 27.9 30.6 33.2 36.3];%摄氏度
5
t=[22.09 16.23 12.65 10.04 8.75];%秒
6
L=[20 20.4 20.4 20.4 20.4]/100;%米
7
D=[8.2 7.88 7.888 7.88 7.88]/100;%米
8
rho_ball=7.700*10^3;
9
rho_oil=0.9741*10^3;
10
g=9.8066;
11
12
13
%计算小球各自直径，用行向量d储存
14
x=x2-x1;%毫米
15
for i=0:4
16
tmp=0;
17
for j=1:5
18
tmp=tmp+x(5*i+j);
19
end
20
d(i+1)=double(tmp/5000);
21
end
22
23
%eta的计算
24
eta=((rho_ball-rho_oil).*g.*(d.^2).*t)./(18.*L.*(1+2.4.*(d./D)));
25
26
%散点图绘制和平滑曲线连接
27
figure(1);
28
T=double(T);
29
T_line=linspace(20,40,2000);
30
eta=double(eta);
31
eta_line=interp1(T,eta,T_line,'cubic');
32
plot(T,eta,"bo",T_line,eta_line,"r-",LineWidth=3);
33
xlabel("温度/℃","FontName","宋体","FontSize",24);
34
ylabel("黏度系数η/Pa·s","FontName","宋体","FontSize",24);
35
title("蓖麻油黏度系数与温度关系图","FontName","宋体","FontSize",24);
36
37
%室温下不确定度的计算
38
syms a b time phi;
39
y=((a-b).*g.*(phi.^2).*time)./(18.*L(1).*(1+2.4.*(phi./D(1))));
40
41
dy=[diff(y,a) diff(y,b) diff(y,time) diff(y,phi)];
42
dy=subs(dy,a,rho_ball);
43
dy=subs(dy,b,rho_oil);
44
dy=subs(dy,time,t(1));
45
dy=subs(dy,phi,mean(d(1:5)));
46
dy=double(dy);
47
48
phi=d(1:5);
49
50
u_a=5;
51
u_b=0.5;
52
u_time=0.2./sqrt(3);
53
u_phi=std(phi,1)*1.14;
54
u_0=[u_a u_b u_time u_phi];
55
u_0=u_0';
56
57
U=sqrt((dy.^2)*(u_0.^2));
58
E=U/eta(1);
59
60
%输出结果
61
disp("小球测量直径依次为(单位：mm）：");
62
x
63
disp("小球平均直径为（单位：mm）：");
64
d.*1000
65
disp("计算得各自黏度（单位：Pa·s）：");
66
eta
67
disp("计算得室温下黏度对小球密度、蓖麻油密度、时间和小球直径各自的一阶偏导为：");
68
dy
69
disp("各自合成不确定度为：");
70
u_0'
71
disp("室温下黏度的总合成不确定度为U="+U+"Pa·s");
72
disp("相对不确定度为E="+E*100+"%");
73
disp("室温下黏度的表达式为η="+eta(1)+"±"+U+"Pa·s");
911%线性和非线性元件伏安特性曲线
2
%国际单位制
3
UR=[1.020 1.200 1.4440 1.704 1.942 2.240 2.900];
4
IR=[0.6137 0.7164 0.8610 1.0211 1.1626 1.3429 1.7401]/1000;
5
UL=[1.5790 1.7448 1.8129 1.8397 1.8656 1.9379 1.9536 2.0714 2.1674 2.2470];
6
IL=[0.0153 0.3267 0.9967 1.3893 1.8805 3.557 4.000 8.075 12.850 19.330]/1000;
7
8
%非线性元件伏安特性曲线绘制
9
figure(1);
10
plot(UL,IL*1000,"b^",LineWidth=2);
11
hold on;
12
U=1.5:0.00005:2.25;
13
I=interp1([1.5 UL],[0 IL],U,'spline');
14
plot(U,I*1000,"r-",LineWidth=3);
15
hold off;
16
xlabel("U/V","FontName","宋体","FontSize",24);
17
ylabel("I/mA","FontName","宋体","FontSize",24);
18
title("发光二极管伏安特性曲线","FontName","宋体","FontSize",24);
19
legend("测量值","拟合值","Location","southeast");
20
grid on;
21
grid minor;
22
23
%线性元件伏安特性曲线绘制
24
figure(2);
25
plot(UR,IR*1000,"b^",LineWidth=2);
26
hold on;
27
m=polyfit(UR,IR,1);
28
u=linspace(0,3.1,2000);
29
i=m(1).*u+m(2);
30
plot(u,i*1000,"r-",LineWidth=3);
31
hold off;
32
xlabel("U/V","FontName","宋体","FontSize",24);
33
ylabel("I/mA","FontName","宋体","FontSize",24);
34
title("恒定电阻伏安特性曲线","FontName","宋体","FontSize",24);
35
legend("测量值","拟合值","Location","southeast");
36
grid on;
37
grid minor;
38
39
%线性元件电阻
40
u_av=0;
41
for n=1:7
42
u_av=u_av+UR(n);
43
end
44
u_av=u_av/7;
45
46
i_av=0;
47
for n=1:7
48
i_av=i_av+IR(n);
49
end
50
i_av=i_av/7;
51
52
u2_av=0;
53
for n=1:7
54
u2_av=u2_av+UR(n).^2;
55
end
56
u2_av=u2_av/7;
57
58
ui_av=0;
59
for n=1:7
60
ui_av=ui_av+UR(n).*IR(n);
61
end
62
ui_av=ui_av/7;
63
64
k=(u_av*i_av-ui_av)/(u_av.^2-u2_av);
65
R=1/k-150;
66
67
%不确定度计算
68
syms z x y;
69
z=x/y;
70
tp=1.09;
71
uu=std(UR,1)*tp;
72
ui=std(IR,1)*tp;
73
%diff=[diff(z,x).^2 diff(z,y).^2];
74
dif=sqrt((1/u_av).^2.*uu.^2+(1/i_av).^2.*ui.^2)*R;
75
what=3*0.005/(sqrt(3)*i_av);%B类不确定度
76
fuck=sqrt(what.^2+dif.^2);%合成不确定度
77
78
%数据输出
79
disp("恒定电阻：");
80
disp("电压平均值："+u_av+"V");
81
disp("电流平均值："+i_av*1000+"mA");
82
disp("电压平方平均值："+u2_av+"V2");
83
disp("电压电流乘积平均值："+ui_av*1000+"mA*V");
84
disp("计算得斜率k="+k*1000);
85
disp("计算得电阻值R="+R+"Ω");
86
disp("计算机斜率"+m(1));
87
disp("计算机截距"+m(2));
88
disp("A类不确定度Ω："+double(dif));
89
disp("B类不确定度Ω："+double(what));
90
disp("合成不确定度Ω："+double(fuck));
91
disp("置信概率：68.3%");
771%示波器的原理和应用
2
%ch1
3
figure(1);
4
T=210;%微秒
5
x=linspace(0,T+30,10000);
6
A=4.4;%V
7
for n=1:10000
8
if(mod(x(n),T)<T/2)
9
y(n)=A;
10
else
11
y(n)=-A;
12
end
13
end
14
plot(x,y,"r-",LineWidth=3);
15
title("CH1:方波",FontName="宋体",FontSize=16);
16
xlabel("t/μs",FontName="宋体",FontSize=14);
17
ylabel("U/V",FontName="宋体",FontSize=14);
18
text(220,-2,['T=210μs' sprintf('\n') 'A=4.4V'],fontname="宋体",FontSize=16);
19
grid on;
20
grid minor;
21
22
%ch2
23
figure(2);
24
T=105;%微秒
25
x=linspace(0,T+30,10000);
26
A=2.3;%V
27
y=abs(A*sin(2*pi/T*x));
28
plot(x,y,"r-",LineWidth=3);
29
title("CH2:全波整流",FontName="宋体",FontSize=16);
30
xlabel("t/μs",FontName="宋体",FontSize=14);
31
ylabel("U/V",FontName="宋体",FontSize=14);
32
text(110,0.5,['T=105μs' sprintf('\n') 'A=2.3V'],fontname="宋体",FontSize=16);
33
grid on;
34
grid minor;
35
36
%ch3
37
figure(3);
38
T=210;%微秒
39
x=linspace(0,T+30,10000);
40
A=4.4;%V
41
y=A*sin((2*pi)/T*x);
42
plot(x,y,"r-",LineWidth=3);
43
title("CH3:正弦波",FontName="宋体",FontSize=16);
44
xlabel("t/μs",FontName="宋体",FontSize=14);
45
ylabel("U/V",FontName="宋体",FontSize=14);
46
text(220,-2,['T=210μs' sprintf('\n') 'A=4.4V'],fontname="宋体",FontSize=16);
47
grid on;
48
grid minor;
49
50
%ch4
51
figure(4);
52
T=210;%微秒
53
x=linspace(0,T+30,10000);
54
A=4.6;%V
55
y=A*sawtooth(2*pi/T*(x+52.5),0.5);
56
plot(x,y,"r-",LineWidth=3);
57
title("CH4:三角波",FontName="宋体",FontSize=16);
58
xlabel("t/μs",FontName="宋体",FontSize=14);
59
ylabel("U/V",FontName="宋体",FontSize=14);
60
text(200,-3,['T=210μs' sprintf('\n') 'A=4.6V'],fontname="宋体",FontSize=16);
61
grid on;
62
grid minor;
63
64
%ch5
65
figure(5);
66
T=210;%微秒
67
x=linspace(0,2*T,10000);
68
A=2.2;%V
69
y=A*sin(2*pi/T*x);
70
y(y<0)=0;
71
plot(x,y,"r-",LineWidth=3);
72
title("CH5:半波整流",FontName="宋体",FontSize=16);
73
xlabel("t/μs",FontName="宋体",FontSize=14);
74
ylabel("U/V",FontName="宋体",FontSize=14);
75
text(350,2,['T=210μs' sprintf('\n') 'A=2.3V'],fontname="宋体",FontSize=16);
76
grid on;
77
grid minor;
551%光的等厚干涉
2
%牛顿环部分
3
4
%计算Dm
5
d2=[27.880 27.415 27.301 27.130 26.989];
6
d1=[16.433 16.535 16.693 16.834 16.997];
7
Dm=(d2-d1)*0.001;%米
8
9
%计算Dk
10
d4=[26.834 26.675 26.536 26.359 26.171];
11
d3=[17.146 17.287 17.446 17.630 17.829];
12
Dk=(d4-d3)*0.001;%米
13
14
%钠光灯波长
15
lambada=589.3.*(10.^-9);%589.3nm
16
17
D2=Dm.^2-Dk.^2;%计算各组数据的平方之差
18
D2ba=mean(D2);%计算平方之差的平均数
19
Rba=D2ba/(20*lambada);%计算平均曲率半径
20
21
%计算不确定度
22
Sd=std(D2,1);%A类不确定度
23
Ud=Sd;%合成不确定度
24
Ur=Ud/(20*lambada);
25
E=Ur/Rba;
26
27
%输出结果
28
disp("曲率半径平均值Rba="+Rba+"m");
29
disp("平方之差不确定度为Ud="+Ud);
30
disp("曲率半径不确定度为Ur="+Ur);
31
disp("曲率半径相对不确定度为E="+E);
32
disp("最终测定，曲率半径表达式为R="+Rba+"±"+Ur+"(m)");
33
disp("各直径依次为（mm）：");
34
D=[Dm Dk]*1000 %毫米
35
36
%空气劈尖部分
37
x0=[2.601 2.605 2.604 2.603 2.602];
38
x1=[35.371 35.373 35.368 35.368 35.370];
39
L=x1-x0;
40
Lba=mean(L);%以上均为毫米
41
42
x2=[8.050 8.661 9.147 9.595 10.064];
43
x3=[12.972 13.443 13.972 14.417 14.920];
44
delta=x3-x2;
45
delta_ba=mean(delta);%以上均为毫米
46
n=10/delta_ba;%一毫米的条纹
47
d=Lba*n*lambada/2;%厚度，米
48
49
%输出结果
50
disp("劈尖长度");
51
L
52
disp("十条宽度");
53
delta
54
disp("一毫米条纹数n="+n);
55
disp("磁带厚度d="+d+"m");