-2

问题说：在铝棒上进行了三个拉伸试验。在每个测试中，应变是在相同的应力值下测量的。结果是 在此处输入图像描述

其中应变单位为 mm/m。使用线性回归估计钢筋的弹性模量（弹性模量 = 应力/应变）。

我用这个程序来解决这个问题：

function coeff = polynFit(xData,yData,m)
% Returns the coefficients of the polynomial
% a(1)*x^(m-1) + a(2)*x^(m-2) + ... + a(m)
% that fits the data points in the least squares sense.
% USAGE: coeff = polynFit(xData,yData,m)
% xData = x-coordinates of data points.
% yData = y-coordinates of data points.

A = zeros(m); b = zeros(m,1); s = zeros(2*m-1,1);
for i = 1:length(xData)
temp = yData(i);
for j = 1:m
b(j) = b(j) + temp;
temp = temp*xData(i);
end
temp = 1;
for j = 1:2*m-1
s(j) = s(j) + temp;
temp = temp*xData(i);
end
end
for i = 1:m
for j = 1:m
A(i,j) = s(i+j-1);
end
end 
% Rearrange coefficients so that coefficient
% of x^(m-1) is first
coeff = flipdim(gaussPiv(A,b),1);

问题在没有程序的情况下解决如下 在此处输入图像描述

我的尝试

T=[34.5,69,103.5,138];

D1=[.46,.95,1.48,1.93];

D2=[.34,1.02,1.51,2.09];

D3=[.73,1.1,1.62,2.12];

Mod1=T./D1;

Mod2=T./D2;

Mod3=T./D3;

xData=T;

yData1=Mod1;

yData2=Mod2;

yData3=Mod3;

coeff1 = polynFit(xData,yData1,2);

coeff2 = polynFit(xData,yData2,2);

coeff3 = polynFit(xData,yData3,2);

x1=(0:.5:190);

y1=coeff1(2)+coeff1(1)*x1;

subplot(1,3,1);

plot(x1,y1,xData,yData1,'o');

y2=coeff2(2)+coeff2(1)*x1;

subplot(1,3,2);

plot(x1,y2,xData,yData2,'o');

y3=coeff3(2)+coeff3(1)*x1;

subplot(1,3,3);

plot(x1,y3,xData,yData3,'o');

我该怎么做才能得到这个结果？

4

1 回答 1

3

As a general advice:

  • avoid for loops wherever possible.
  • avoid using i and j as variable names, as they are Matlab built-in names for the imaginary unit (I really hope that disappears in a future release...)

Due to m being an interpreted language, for-loops can be very slow compared to their compiled alternatives. Matlab is named MATtrix LABoratory, meaning it is highly optimized for matrix/array operations. Usually, when there is an operation that cannot be done without a loop, Matlab has a built-in function for it that runs way way faster than a for-loop in Matlab ever will. For example: computing the mean of elements in an array: mean(x). The sum of all elements in an array: sum(x). The standard deviation of elements in an array: std(x). etc. Matlab's power comes from these built-in functions.

So, your problem. You have a linear regression problem. The easiest way in Matlab to solve this problem is this: 

%# your data
stress = [ %# in Pa
    34.5 69 103.5 138] * 1e6;

strain = [ %# in m/m
    0.46 0.95 1.48 1.93
    0.34 1.02 1.51 2.09
    0.73 1.10 1.62 2.12]' * 1e-3;    

%# make linear array for the data
yy = strain(:);
xx = repmat(stress(:), size(strain,2),1);

%# re-formulate the problem into linear system Ax = b
A = [xx ones(size(xx))];
b = yy;

%# solve the linear system
x = A\b;

%# modulus of elasticity is coefficient 
%# NOTE: y-offset is relatively small and can be ignored)
E = 1/x(1)

What you did in the function polynFit is done by A\b, but the \-operator is capable of doing it way faster, way more robust and way more flexible than what you tried to do yourself. I'm not saying you shouldn't try to make these thing yourself (please keep on doing that, you learn a lot from it!), I'm saying that for the "real" results, always use the \-operator (and check your own results against it as well).

The backslash operator (type help \ on the command prompt) is extremely useful in many situations, and I advise you learn it and learn it well.

I leave you with this: here's how I would write your polynFit function: 

function coeff = polynFit(X,Y,m)

    if numel(X) ~= numel(X)
        error('polynFit:size_mismathc',...  
              'number of elements in matrices X and Y must be equal.');
    end

    %# bad condition number, rank errors, etc. taken care of by \
    coeff = bsxfun(@power, X(:), m:-1:0) \ Y(:);

end

I leave it up to you to figure out how this works.

于 2012-10-08T06:48:53.907 回答