Mathstudio在大学物理实验数据处理中应用
2024-06-20
摘要:将数学软件Mathstudio应用到大学物理实验数据处理中,进行描述性统计、推断统计、求不确定度、线性回归等运算。Mathstudio具备数值运算和符号运算功能,使用数组和切片(Slice)操作,内置大量数学函数,微积分、统计等功能很强大,作图和动画也方便。Mathstudio不用安装、编译,浏览器打开网址即可运行,可逐行调试,命令格式简单。示例结合线性代数理论,使用了雅可比矩阵、海森矩阵、范数、线性回归、作图等命令,实现Mathstudio编程计算空心圆柱体体积的不确定度、铜-康铜热电偶温差电势的线性回归模型,程序简短精练,结构清晰,提高了数据处理效率。Mathstudio编程效率高,难度较低,适合小规模数据快速分析,也能进一步开发更专业的数据处理功能。
关键词:Mathstudio 描述统计 推断统计 梯度 不确定度 线性回归
中图分类号:G633.7 文献标识码:A
Application of Mathstudio in the Experimental Data Processing of University Physics
ZHOU Hongliang
(Jiangsu Vocational College of Electronics and Information, Huai'an, Jiangsu Province, 223003 China)
Abstract: Mathematical software Mathstudio is applied to the experimental data processing of university physics to perform operations such as descriptive statistics, inferential statistics, uncertainty and linear regression. Mathstudio has the function of numerical and symbolic operations, uses arrays and slice operations, has a large number of built-in mathematical functions, has powerful functions such as calculus and statistics, and is also convenient for graphing and animating. Mathstudio does not need to be installed and compiled, and it can be run by opening the Web site in the browser and debugged line by line with the simple command format. Combined with the linear algebra theory, the example uses commands such as the Jacobi matrix, Hessian matrix, paradigm, linear regression and graphing to realize the calculation of the uncertainty of the volume of hollow cylinders and the linear regression of copper-constantan thermocouple temperature difference potential by Mathstudio programming, and the program is short and concise with clear structure, which improves the efficiency of data processing. Mathstudio programming is highly efficient and less difficult, and it is suitable for the rapid analysis of small-scale data and also can further develop the more professional functions of data processing.
Key Words: Mathstudio; Descriptive statistics; Inferential statistics; Gradient; Uncertainty; Linear regression
大学物理实验数据处理与统计密切相关,针对不同实验目的,实验数据处理包括描述性统计、不确定度分析、线性拟合、回归分析、向量微分、导数、梯度等运算。教材中的作图法、逐差法、最小二乘法等等是实验数据处理的理论基础,若采用手工计算,过程繁琐并且错误率高,使用数学软件可以有效提高数据处理效率,便于实现数据可视化。
Mathstudio是一款小巧精致的数学软件[1],无须编译,在网页http://mathstud.io/输入指令直接运行,可以进行统计运算、微积分运算、级数运算、傅里叶变换、矩阵运算(特征值、QR分解、LU分解、SVD分解)、解方程组和微分方程等,函数与MATLAB类似,有丰富的作图和动画功能,提供循环判断语句和逻辑关系等运算,可以实现较完整的功能。在矩阵运算上相当便捷,可以使用函数直接对矩阵计算,命令简洁,比循环语句效率高。
1描述性统计
1.1集中量数
1.1.1向量计算
设随机变量,X容量为n,常用的统计集中量有计数、求和、最大值、最小值、平均值、中位数、众数。
Mathstudio命令[length(X), sum(X), max(X), min(X), mean(X), median(X)]计算各统计量(无mode函数),与MATLAB命令格式相同。
1.1.2矩阵计算
设随机变量有m组样本,用矩阵表达,第i列向量。示例:Mathstudio输入矩阵X = [[5, 0, 3, 7], [1, -5, 7, 3], [4, 9, 8, 10]],取第1列命令X(*, 1),得到列向量。
Mathstudio命令[max(X), min(X)]统计矩阵所有元素的最值,[mean(X), median(X)]按矩阵列方向计算并返回行向量。与Matlab区别在于,Matlab除length()以外的函数都是按矩阵列计算返回行向量。
1.2差异量数
1.2.1标准差和方差
求标准差(Standard Deviation)和方差(Variance)是数据分析的基础,众多统计量数包含标准差或形式上类似标准差的计算。容量为n的随机变量X,X的标准差σ用于度量个体xi和X的数学期望E(X)(即X的均值)之间的偏离程度,或者说表达了随机变量数据的离散程度,方差是标准差的平方。
样本标准差分无偏估计和有偏估计,分母分别取n?1(当n>1时)和n计算。
Mathstudio命令StandardDeviation(X)计算向量X的无偏估计标准差,或者对矩阵X列向量计算并返回行向量。Variance(X)计算无偏估计标准方差。若要按行计算用StandardDeviation(X'),其中撇号“'”表示转置。
1.2.2标准误差和不确定度
样本均值的标准误差(Standard Error),,描述样本均值的分布,在物理实验中用于表达测量列X的A类不确定度,测量列X的A类不确定度uA = StandardDeviation(X)/SQRT(Length(X))。
2推断统计
2.1协方差
协方差(Covariance)用来衡量随机变量X、Y的总体误差,表达了X和Y的相关性,
。矩阵X的协方差矩阵
是对称、半正定矩阵,对角线为两个变量相同时的协方差,即。
2.2相关系数
皮尔逊相关系数[2](Pearson correlation coefficient),反映线性相关变量之间相关密切程度的统计指标,相关系数将协方差归一化(normalization)消除变量量纲或者变化幅度不同的影响,单纯反映两个变量在每单位变化上的相似程度,,其中,,余类似。
,越接近1说明X、Y线性相关越明显,如果接近0说明X、Y不相关,不能用线性回归模型来估计X、Y的关系。
2.3范数
向量X的2范数,即取向量的模,在三维坐标系中2范数表示点M(x, y, z)到原点O的欧几里得距离(Euclidean distance)。在物理实验中求间接测量不确定度时,将各不确定度分量合成,需要求元素的平方和再开根号的运算,就可使用norm()。
评定测量结果的最小二乘估计模型的精度时,表达模型误差通常用残差平方和,即残差向量ε的2范数的平方,或者使用均方根误差。
Mathstudio中Norm(X)返回向量的2范数,模型误差可用RSS = Norm(ε)^2和RMSE = Norm(ε) / SQRT(Length(ε))计算。
2.4 线性回归
某实验测量值(xi, yi)分属于和,若物理量y和x之间存在线性关系,线性回归(Linear Regression)指找到两变量y、x之间的线性相关性,有时也称线性拟合(Linear Fit)。根据最小二乘法[3]的原理,找到一条最佳的拟合曲线,那么各测量值与这条拟合曲线上对应点之差(即残差)的平方和为最小。设估计的一元线性回归方程,其中w为斜率,b为截距。设定目标函数并要求J最小,即J→min,可见J是w和b的函数。求J对w和b的偏导数并令偏导数等于0,解得,,就是线性回归方程中的待定参数w和b的最佳估计值。
在待定参数w和b确定以后,为了判断所得结果是否合理,通常用相关系数r来检验[4],,值越接近1,说明实验数据点密集地分布在所拟合的直线附近,用最小二乘法进行线性回归是合适的。时可认为两个物理量之间存在较密切的线性关系,线性回归才有意义。表示变量X、Y完全线性相关,拟合直线通过全部实验数据点。相反,值越小线性越差,如果而接近于0,用线性回归来推测实验数据的模型与实际差异很大,说明假设是错误的,应采用其他函数曲线或方法进行拟合。
Mathstudio用LinearFitModel()计算线性回归。
3 导数和微分
在计算不确定度、线性拟合时需要运用一阶偏导数和二阶偏导数运算。
3.1雅可比矩阵
对列向量,有行向量偏导算子,其作用于函数向量,得到雅可比(Jacobian)矩阵J,函数按列方向,自变量按行方向,
3.2梯度矩阵
对列向量,有列向量偏导算子,其作用于函数向量,得到梯度(Gradient)矩阵:
梯度矩阵是其雅可比矩阵的转置。梯度与方向导数关系密切,方向导数的最大值为梯度向量的模。在梯度正方向,函数以最大速率上升。例如电场强度的方向是负梯度方向,沿着电场线方向电势下降最快。
3.3 海森矩阵
如果实值多元函数二阶连续可导,海森(Hessian)矩阵是梯度对自变量x的Jacobian矩阵,
在判断多元函数极值时,在驻点M有,是可能的极值点。但是仅通过一阶偏导数无法判断,需要通过二阶偏导数判断极值。
Hessian矩阵为二阶偏导数矩阵,为实对称阵,根据其对应的二次型,Hessian矩阵有正定、半正定、负定、半负定、不定等类型。Hessian矩阵正定性与多元函数的凹凸性有关。
正定[5],M为极小值点;负定,M为极大值点;不定矩阵,M非极值点;为半正定矩阵(或半负定矩阵)时,M为“可疑”极值点,需要利用其他方法来判定。
Mathstudio分别用Jacobian()、Gradient()、Hessian()计算雅可比矩阵、梯度矩阵、海森矩阵。
4 Mathstudio实验数据处理示例
4.1 计算空心圆柱体积不确定度
4.1.1 数据和变量预处理
clear(all) // 清空变量
vars = [D, d, H, h] // 定义符号自变量外径D、内径d、高度H、深度h,单位mm
V = 3.1416 / 4 * (D^2 * H - d^2 * h) // 定义体积符号函数
D0 = [16.00, 16.01, 16.02, 16.05, 15.93, 15.99]
d0 = [14.02, 14.03, 14.02, 14.08, 14.22, 14.01]
H0 = [20.01, 20.02, 20.00, 20.06, 20.10, 19.88]
h0 = [18.01, 18.22, 18.00, 18.06, 18.14, 18.18]
data = [D0, d0, H0, h0]' // 原始数据矩阵转置
zero_offset = [0, 0, 0, 0] // 长度测量零点误差
mean_list = mean(data) - zero_offset // 外径、内径、高度、深度平均值列表
uA_list = StandardDeviation(data) / SQRT(length(data(*, 1))) // A类不确定度列表
uB = [0.02, 0.02, 0.02, 0.02] / SQRT(3) // 游标卡尺B类不确定度
uC_list = sqrt(uA_list^2 + uB^2) // C类不确定度列表
说明:程序用数组Data存储测量原始数据,对外径D、内径d、高度H和深度h定义符号自变量vars = [D, d, H, h]和符号函数,第二步Mean(X)计算自变量的平均值,并用标准误差计算A类不确定度,计算B类不确定度,再用计算合成不确定度,全部以矩阵函数f(data)直接计算。
4.1.2计算体积结果和间接测量不确定度
V_val = eval(V, vars, mean_list) // 体积结果
Jaco_mat = Jacobian(V, vars) // 求体积偏导函数
Jaco_vals = Jacobian(V, vars, mean_list) // 求体积偏导函数并赋值
uV = Norm(Jaco_vals * uC_list) // 总不确定度
说明:经间接测量得各分量的不确定度,体积的间接测量不确定度[6]
,如果手工计算非常繁琐,极易出错。将上式改写成,运用符号运算Jacobian()计算体积的雅可比矩阵,求出梯度向量gradV,gradV和取哈达玛积(Hadamard product)得到偏微分向量,再用norm()取向量的模,即平方和开根号,计算体积的总不确定度。
程序在均值、标准差、范数等运算中,将变量定义为向量或矩阵,函数可以作用于矩阵A,直接用矩阵函数,命令含义一目了然。矩阵运算功能简化了程序,避免使用循环,程序结构更清晰,运算效率提高。利用Mathstudio强大的符号运算功能计算偏微分,精简了语句。
4.1.3 数据显示
m = 10^floor(log(10, uV)) // uV保留一位缩放
res = ([round(V_val / m), ceil(uV / m)] * m) // 体积四舍五入,不确定度只进不舍
n = floor(log(10, V_val)) // 科学计数法n次幂
res = res / 10^n // 科学计数法,除以10^n
["V=", res(1), "±", res(2), "×10^", n] // 显示结果V_val±uV=1.21±0.02×10^3
说明:10^floor(log(10, uV))计算出采用科学计数法需要缩放的倍率,不确定度的保留规则是只进不舍,保留一位,与体积V有效数字右对齐,uV用ceil()进位,V用round()四舍五入,结果用科学计数法显示。
4.1.4 梯度矩阵、海森矩阵极值判定
fDiff(V, vars) // 计算全微分,求偏导函数
grad_V = Gradient(V, vars) // 求梯度函数
grad_vals = eval(grad_V, vars, mean_list) // 梯度函数赋值
stationary_point = SolveSystem(grad_V, vars) // 计算一阶偏导数为0的点
Hess_mat = Hessian(V, vars) // 计算海森矩阵
Hess_vals = eval(Hess_mat, vars, stationary_point(2,*)) // 符号变量赋驻点的值
eig_vals = Eigenvalues(Hess_vals) // 计算驻点海森矩阵特征值,判断正负定
说明:对比fDiff()、Gradient()、Jacobian()三个函数,都可以计算偏导数,fDIff()是全微分,可以简记全微分为,Gradient()和Jacobian()的结果在Mathstudio中都用偏导数向量表示。
在求出偏导数后,根据多元函数极值的条件,求解梯度方程组,SolveSystem(grad_V, vars)计算V对[D, d, H, h]的一阶偏导数为0的点(注:结果仅供参考),Hessian()计算海森矩阵即V对[D, d, H, h]的二阶偏导数矩阵,eval()将驻点自变量值赋值给海森矩阵,Eigenvalues()计算驻点处海森矩阵的特征值,由特征值判断正、负定,判断极值点。本例为不定型,无法判断极值,需要其他手段再判断。
4.2 热电偶温差系数线性回归
4.2.1 温差电动势
热电偶的温差电动势大小由热端和冷端的温差决定,其热端为正极,冷端为负极,当已知冷端温度,并测出其温差电动势后,便可求出热端温度。温差电动势的大小与金属材料的性质及接触点处的温度差有关。当冷热端温差不大时,温差电动势与温度的关系近似为,式中:ε为温差电动势;T1为热端温度;T0为冷端温度,a是由构成热电偶的金属材料决定的常数。
温差电动势ε与冷热端温差ΔT=T1?T0成线性关系,比率系数a称为热电偶的温差系数或称热电动势率。而在一般情况下,温差电动势与两接触点温差的关系为一曲线。
铜-康铜热电偶的温差系数近似地取a=0.0435 mV/K。但由于组成热电偶的材料成分有差异,所以实际数值有出入。因此在实际使用时,要先对热电偶进行校正,以确定a的数值。校正时,热端缓慢加热,并测量热端温度T1,冷端保持固定T0,可以在杜瓦瓶中放冰水混合物,测出不同温度差下的温差电动势,将ε和ΔT运用线性回归,求出斜率和截距,即得到温差系数。
4.2.2铜-康铜热电偶的温差系数线性回归
T1 = [30, 40, 50, 60, 70, 80, 90, 100, 110, 120] // 热电偶热端温度,单位℃
T0 = 19 // 热电偶冷端温度,单位℃
dT = T1 - T0 // 温度差ΔT
ε = [0.53, 0.94, 1.36, 1.81, 2.25, 2.70, 3.15, 3.62, 4.08, 4.56] // 温差电势ε,单位mV
m = LinearFitmodel(dT, ε) // 线性回归模型斜率m.a、截距m.b和相关系数m.r
说明:程序传递参数ΔT、ε两个向量作为线性拟合的参数,利用LinearFitmodel()函数直接返回线性回归模型,包括最佳估计的斜率w、截距b、相关系数r三个参数,求得温差电动势与温度差的回归模型,温度系数0.0449mV/K,相关系数。
4.2.3 回归结果分析
Linearfitplot(dT, ε) // 线性回归作图
err = m.a * dT + m.b - ε // 残差向量
RMSE = SQRT(Mean(err^2)) // 模型均方根误差
说明:相关系数接近1,表明ε与ΔT强线性相关,回归模型使用得当。Linearfitplot()可以作出数据的散点图和拟合直线对比,从图2中直观看出数据点与拟合直线非常接近,计算结果残差的均方根误差RMSE=0.026,印证了回归模型的精度。温差系数的计算结果0.044 9 mV/K与理论值0.043 5 mV/K非常接近,相对误差3%,测量精度高,数据处理效果很好。
5 结论
Mathstudio小巧灵活,容易入门,提供大量的内置数学函数包括统计、导数微分、符号运算、矩阵运算和作图等,具备关系、逻辑运算和循环、判断语句,运用数组切片操作和矩阵函数,语句简洁,程序结构清晰,避免了循环判断语句,提高数据处理效率。运用Mathstudio编程进行大学物理实验不确定度和线性回归计算,具有方便快捷的优点,易于修改应用到其他实验数据的分析处理。Mathstudio适合大学物理实验数据的处理,对促进教学数字化转型,提升学生自主学习能力、信息素养与持续发展能力很有意义。
参考文献
[1]蔡春雨,萨仁高娃,包琳,等.Mathstudio在近代物理实验数据处理中的应用[J].大学物理实验,2023,36(2):139-143.
[2]李阳,李晓琪,杨佳滢,等.基于THz时域反射成像技术的玉米种子淀粉分布可视化研究[J].光谱学与光谱分析,2023,43(9):2722-2728.
[3]吴城,李晓娟.基于单线激光的果园机器人自主导航方法[J].江苏大学学报(自然科学版),2023,44(5):530-539.
[4]曹久莹,于陆军.基于最小二乘法拟合的流量计不确定度分析方法[J].中国测试,2022,48(S1):122-128.
[5]田宇洋,顾青涛,张任楠,等.协同定位最小二乘凸性分析[J].现代电子技术,2022,45(3):10-16.
[6]贾欣.原子吸收分光光度计石墨炉法测定小米中铬的不确定度分析[J].食品工程,2023(3):44-47.