用matlab实现误差理论与数据处理的部分习题

误差理论与数据处理作业第一章Matlab代码：%题1-5%180度2分换算为clear all;clc;test5 = 180 + 2/60/60 ;absoluteError5 = (180 - test5)oppositeError5 = absoluteError5 / test5 * 100 %%1-7absouluteError7 = 100.2 -100.5oppositeError7 = absouluteError7 / 100.5 * 100%题目1-14比较测量方式的精度L1 = 110;L2 = 110;testError1 = 0.011;testError2 = 0.009;L3 = 150;testError3 = 0.012;%相对误差越小则越精确oppositeError1 = testError1 / L1;oppositeError2 = testError2 / L2;oppositeError3 = testError3 / L3;%寻找相对误差的最小值if oppositeError1 > oppositeError2if oppositeError2 > oppositeError3fprintf('精度最高的是第三种，相对误差为 %f',oppositeError3); elsefprintf('精度最高的是第二种，相对误差为 %f',oppositeError2); endelseif oppositeError1 > oppositeError3fprintf('精度最高的是第三种，相对误差为 %f',oppositeError3); elsefprintf('精度最高的是第一种，相对误差为 %f',oppositeError1); endend运行结果：absoluteError5 =-5.5556e-04oppositeError5 =-3.0864e-04absouluteError7 =-0.3000oppositeError7 = -0.2985精度最高的是第三种，相对误差为 0.000080第二章Matlab代码：这章自己写了几个函数：function value = standardDeviation(n,v)%value是返回的标准差,n 是测量的值的矩阵，v是算术平均值% standardDeviation函数是求测量的标准差k = 1;for i = nresidualError(k,:) = i - v ; %存放残差的向量k = k+1;endvalue = ((sum(residualError.^2)) / (numel(n) - 1))^0.5;endfunction value = T_fenbubiao(a,v)%该函数实现T分布表的查询%v表示次数，a表示可信度m_005 = [12.71 4.30 3.18 2.78 2.57 2.45 2.36 2.31 2.26 2.23;2.20 2.18 2.16 2.14 2.13 2.12 2.11 2.10 2.09 2.09;2.08 2.07 2.07 2.06 2.06 2.06 2.05 2.05 2.05 2.04;0 0 0 2.02 2.01 2.00 1.99 1.99 1.99 1.98 ];m_001 = [63.66 9.92 5.84 4.60 4.03 3.71 3.50 3.36 3.25 3.17;3.11 3.05 3.01 2.98 2.95 2.92 2.90 2.88 2.86 2.85;2.83 2.82 2.81 2.80 2.79 2.78 2.77 2.76 2.76 2.75;0 0 0 2.70 2.68 2.66 2.65 2.64 2.63 2.63 ];m_00027= [235.80 19.21 9.21 6.62 5.51 4.90 4.53 4.28 4.09 3.96;3.85 3.76 3.69 3.64 3.59 3.54 3.51 3.48 3.45 3.42;3.40 3.38 3.36 3.34 3.33 3.32 3.30 3.29 3.28 3.27;0 0 0 3.20 3.16 3.13 3.11 3.10 3.09 3.08 ]; p = a ;switch pcase double(0.05)index = lookup_index(v);row = index(1);column = index(2);if (row > 0 && column > 0)value = m_005(row,column);elsevalue = 1.96;endcase 0.0100;index = lookup_index(v);row = index(1);column = index(2);if (row > 0 && column > 0)value = m_001(row,column);elsevalue = 2.58;endcase 0.0027;index = lookup_index(v);row = index(1);column = index(2);if (row > 0 && column > 0)value = m_00027(row,column);elsevalue =3.00;endotherwisevalue = 1.96;endendfunction value= lookup_index(v)if v <= 10value = [1,v];elseif (v > 10 && v <= 20)v1 = rem(v,10);if v1 == 0;v1 = 10;endvalue =[2,v1];elseif (v > 20 && v <=30 )v2 = rem(v,10);if (v2 == 0)v2 = 10;endvalue =[3,v2];elseif(v > 30 && v <= 100)v3 = fix(v/10);value = [4,v3];elsevalue = [0,0];endendfunction residualError = calculationResiduals(testValue) %求残差的函数% testValue是传入的测量值，residualError返回的残差值k = 1;averageValue = sum(testValue) / numel(testValue);for i = testValuetemp = i - averageValue;residualError(:,k) = temp ; %存放残差的向量k = k+1;endendfunction value = judge( v,d3 )%判断有没有粗大误差函数% d3是3倍的标准差temp = abs(v);vMax = max(temp);if vMax > d3value = findIndex(temp,vMax);elsevalue =false;endendfunction value = findIndex(v,data)k = 1;for i = vif i == data;value = k;return ;endk = k+1;endvalue = false;end实现题目的主要代码：c lear all;clc;%题目2-4 求平均误差和算数平均值test = [236.45,236.37,236.51,236.34,236.39,236.48,236.47,236.40 ];%8次测试值arithmeticAverageValue = sum(test) / numel(test) ; %算术平均值%根据贝塞公式求单次测量的标准差,d 为标准差，dd为算术平均标准差d = standardDeviation(test,arithmeticAverageValue);dd = d / (numel(test))^0.5;fprintf('题目2-4 ');fprintf('算术平均值为 %f ，单次测量的标准差为 %f,算术平均标准差为 %f \n',arithmeticAverageValue,d,dd);%题目2-6 求算术平均值，标准差，或然误差，平均误差testValue2_6 = [168.41,168.54,168.59,168.40,168.50]; arithmeticAverageValue2_6 = sum(testValue2_6) / numel(testValue2_6); %算术平均值d2_6 = standardDeviation(testValue2_6,arithmeticAverageValue2_6); %标准差dd2_6 = d2_6 / (numel(testValue2_6))^0.5;p2_6 = 0.6745 * dd2_6; %或然误差averageError2_6 = 0.7979 * dd2_6; %平均误差fprintf('题目2-6 ');fprintf('算术平均值为 %f ，单次测量的标准差为 %f,算术平均标准差为 %f \n ',arithmeticAverageValue2_6,d2_6,dd2_6);fprintf(' 或然误差 %f，平均误差 %f \n',p2_6,averageError2_6);%题目2-8 求置信限averageDeviation = 0.005 / (5^0.5);n = 5;P = 0.95; %置信限95%% （1）按t分布计算a = 0.05;v = 4;t = T_fenbubiao(a,v); %根据t分布表查对应的t值x = t * averageDeviation;fprintf('题目2-12');fprintf('t分布下置信限为 -%f 到 %f\n ',x,x);% (2)按正太分布做；t = 1.96;%根据表知t = 1.96x = t * averageDeviation;fprintf('正太分布下置信限为 -%f 到 %f \n',x,x);%题目2-12 求加权算术平均值及其标准差xp = [102523.85 102391.30 102257.97 102124.65 101991.33 101858.01 101724.69 101591.36]; %测试的次数p = [1 3 5 7 8 6 4 2 ];averageXp = sum(xp) / numel(xp);k = 1;for i = xptemp = i - averageXp;residualError(:,k) = temp ; %存放残差的向量k = k+1;endresidualErrorV2 = residualError.^2;sum2_12 = sum(xp.*p);weightedAverage = sum2_12 / sum(p); %加权算术平均值value = (sum(p.*residualErrorV2)/((numel(xp)-1)*sum(p))).^0.5;%加权算术标准差fprintf('题目2-12');fprintf('加权算术平均值 %f ，加权算术标准差 %f\n',weightedAverage,value);%题目2-17 判断是否有系统误差testValue2_17 = [14.7 15 15.2 14.8 15.5 14.6 14.9 14.8 15.1 15]; averageValue = sum(testValue2_17) / numel(testValue2_17);%判断测量次数的奇偶cnt = numel(testValue2_17);if mod(cnt,2) == 0K = cnt/2;elseK = (cnt + 1) / 2;endv = calculationResiduals(testValue2_17);d2_17 = sum(v(1:K)) - sum(v((K+1):cnt));fprintf('题目2-17');if abs(d2_17) > abs(max(v))fprintf('%f 显著不为零,故测量中含有线性误差 ',d2_17);elsefprintf('%f 不显著为零,故测量中不含有线性误差 ',d2_17);end%2-22莱以特准则判断粗大误差testValue2_22 = [28.53 28.52 28.50 28.52 28.53 28.53 28.50 28.49 28.49 28.51 28.53 28.52 28.49 28.40 28.50];n = numel(testValue2_22);averageValue2_22 = sum(testValue2_22) / n ;d2_22 = standardDeviation(testValue2_22,averageValue2_22) ; %计算标准差v2_22 = calculationResiduals(testValue2_22); %残差d3_2_22 = d2_22 * 3 ;while 1index = judge(v2_22,d3_2_22);if index %如果有粗大误差剔除fprintf('\n值%f是粗大误差应该提除\n',testValue2_22(index)); testValue2_22(index) = [];%删除对应的值%重复第一次的步骤n = numel(testValue2_22);averageValue2_22 = sum(testValue2_22) / n ;d2_22 = standardDeviation(testValue2_22,averageValue2_22) ; %计算标准差v2_22 = calculationResiduals(testValue2_22); %残差d3_2_22 = d2_22 * 3 ;elsebreak;endend代码的运行结果：题目2-4 算术平均值为 236.426250 ，单次测量的标准差为 0.059746,算术平均标准差为 0.021124题目2-6 算术平均值为 168.488000 ，单次测量的标准差为 0.082280,算术平均标准差为 0.036797 或然误差 0.024819，平均误差 0.029360题目2-8 t分布下置信限为 -0.006216 到 0.006216 正太分布下置信限为-0.004383 到 0.004383题目2-12 加权算术平均值 102028.342500 ，加权算术标准差 87.661818题目2-17 0.800000 显著不为零,故测量中含有线性误差题目2-22 值28.400000是粗大误差应该踢除第三章MATLAB代码：clear all;clc;%题目3-1 求测量误差test = [40 12 1.25 1.005];testOffset = [-0.7 0.5 -0.3 0.1];limL = [0.35 0.25 0.20 0.20];testSum = sum(test);systemOffset = sum(testOffset);offsetValue = testSum - systemOffset * 1e-3;limLOffset = sum(limL.^2)^0.5;fprintf('题目3-1 ');fprintf('系统误差 %f 修正值为 %f 测量误差 %f \n',systemOffset,offsetValue,limLOffset);%题目3-2 求立方体体积及其极限误差a = 161.6;b = 44.5;c = 11.2;da = 1.2;db = -0.8;dc = 0.5;%都是正负多少lima = 0.8;limb = 0.5;limc = 0.5;V = a*b*c;Vd = (da*b*c) + (a*db*c) + (a*b*dc);Vr = V - Vd;limV = ((b*c)^2*lima^2 + (a*c)^2*limb^2 + (a*b)^2*limc^2).^0.5; fprintf('题目3-2 测量体积最后结果表示为 %d +- %d \n',Vr,limV);%题目3-4 求功耗global U;global I;%global f;syms U IP = U*I;fu = diff(P,U); %求偏导fi = diff(P,I);I = 22.5;U = 12.6;di = 0.5;du = 0.1;ss = vpa(subs(fu));zz = vpa(subs(fi));Dp = (ss^2 * du^2 + zz^2 * di^2 + 2*ss*zz*di*du)^0.5;ui = vpa(subs(P));fprintf('题目3-4 功率P为 %sym 标准差为 %sym \n',ui,Dp);%题目3-12 求误差global r;global h;syms r h%相对误差为1%,测量项目为2项 n = 2V3_12 = pi * h * r^2;Vr = diff(V3_12 ,r);Vh = diff(V3_12,h);n = 2^0.5;e = 0.01;r = 2.0;h = 20.0;V0 = vpa(subs(V3_12));dv = V0 * e;rr = vpa(subs(Vr));hh = vpa(subs(Vh));dr = (dv / (n))*(1 / rr);dh = (dv / (n))*(1 / hh);ddr = vpa(dr);ddh = vpa(dh);fprintf('题目3-12 测定r的误差应为 %0.6sym 测定h的误差应为 %0.6sym \n',vpa(ddr,4),vpa(ddh,5));代码的运行结果如下：题目3-1 系统误差 -0.400000 修正值为 54.255400 测量误差 0.514782题目3-2 测量体积最后结果表示为 7.779570e+04 +- 3.729111e+03题目3-4 功率P为 283.5ym 标准差为 8.55ym题目3-12 测定r的误差应为 0.0070 测定h的误差应为 0.1414第四章MATLAB代码：clear all;clc ;%题目 4-1 求圆球的最大截面的圆周和面积及圆球的体积的测量不确定度%（1）求圆球的最大截面的圆周的测量不确定度global r;syms r;%求导D = 2*pi*r;f = diff(D,r);S = pi * r^2 ;f1 = diff(S,r);V = (4/3)*pi*r^3;f2 = diff(V,r);fprintf('题目 4-1 \n')dr = 0.005;%标准差r = 3.132;fv = vpa(subs(f));Uc = (fv^2*dr^2)^0.5; %标准不确定度n = 9;%测量次数减一t = 0.01;%置信度K = T_fenbubiao(t,n);U = K * Uc; %最大截面的圆周的测量不确定度fprintf('最大截面的圆周的测量不确定度为 %0.5s cm\n',U);%(2) 求圆球的最大截面的面积的测量不确定度f1v = vpa(subs(f1));Uc1 = (f1v^2 * dr^2)^0.5;U1 = K*Uc1;fprintf('最大截面的圆面积的测量不确定度为 %0.6s cm\n',U1); %(3)圆球的体积的测量不确定度f2v = vpa(subs(f2));Uc2 = (f2v^2 * dr^2)^0.5;U2 = K*Uc2;fprintf('圆球的体积的测量不确定度为 %0.6s cm\n',U2);% 题目4-2 求不确定度分量和放大率D的标准不确定度global y1 ;global y2 ;syms y1 y2D = y1 / y2;Uy1 = 0.10;Uy2 = 0.005;Udy1 = diff(D,y1);Udy2 = diff(D,y2);y1 = 19.80;y2 = 0.80;Udy1v = vpa(subs(Udy1))^2;Udy2v = vpa(subs(Udy2))^2;Ud = (Udy1v * Uy1^2 + Udy2v * Uy2^2)^0.5;fprintf('题目4-2 \n');value = vpa(subs(D));fprintf('望远镜放大率D的合成标准不确定度为 %0.6s \n',Ud); fprintf('望远镜放大率D为 %0.6s +- %0.6s (cm)\n ',value,Ud);% 题目4-5 求量块组引起的测量不确定度testValue4_5 = [40 10 2.5];dt = [0.45 0.30 0.25];P = 99.73;ta = 3;%根据P的值查表得Ul = dt./3 ;U = (sum(Ul.^2)).^0.5;data = roundn(U,-3);%保留3位小数fprintf('题目4-5 \n');fprintf(' L的测量不确定度为 %f um \n',data);代码运行结果：题目 4-1最大截面的圆周的测量不确定度为 0.102 cm最大截面的圆面积的测量不确定度为 0.3197 cm圆球的体积的测量不确定度为 2.0031 cm题目4-2望远镜放大率D的合成标准不确定度为 0.1988望远镜放大率D为 24.75 +- 0.1988 (cm)题目4-5L的测量不确定度为 0.1990 um11 / 11。