2012-10-08 17 views
-2

問題所在: 在鋁棒上進行三次拉伸試驗。在每次測試中,應變都是在相同的應力值下測量的。結果爲 enter image description here在matlab中的拉伸試驗

其中應變單位爲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); 

如下 enter image description here

我嘗試的問題是沒有一個程序解決

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'); 

我有什麼做的得到這個結果?

+0

你能更具體地說明問題出在哪裏嗎?你能否澄清'沒有程序的工作'解決方案如何與你正在嘗試的多項式曲線擬合方法相關?你究竟在尋找什麼結果?我在示例中看到很多mean(x),但它在代碼中根本沒有使用......? – Dan

回答

3

作爲一般的建議是:

  • 避免循環儘可能。
  • 忌用ij作爲變量名,因爲它們是Matlab的內置的虛數單位名稱(我真的希望在未來的版本中會消失......)

由於m被一種解釋語言,for循環可能會比編譯後的代碼慢很多。 Matlab被命名爲MATtrix LABoratory,這意味着它針對矩陣/陣列操作進行了高度優化。通常情況下,當有一個操作無法在沒有循環的情況下完成時,Matlab有一個內置函數,它的運行方式比Matlab中的for-loop更快。例如:計算數組中元素的平均值:mean(x)。數組中所有元素的總和:sum(x)。數組中元素的標準偏差:std(x)。 Matlab的功能來自這些內置功能。

所以,你的問題。你有一個線性迴歸問題。在Matlab解決這個問題最簡單的方法是這樣的:

%# 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) 

你在函數polynFit沒有什麼是A\b完成,但\ - 運算符是能夠做到方式更快,方式更多健壯和方式比你試圖做自己更靈活。我並不是說你不應該嘗試自己製造這些東西(請繼續這樣做,你會從中學到很多東西!),我說的是對於「真實」的結果,請始終使用\-operator(並檢查你自己的結果)。

反斜槓操作符(在命令提示符下鍵入help \)在許多情況下非常有用,我建議您學習並學習它。

我離開你這個:這裏是我會怎麼寫你polynFit功能:

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 

我把它留給你找出如何工作的。