%《MATLAB及其在理工課程中的應用指南》(十一五規劃版)第七章例7－2－3程序
% 簡支梁的變形計算
% 西安電子科技大學出版社出版  陳懷琛編著 2007年5月
%
%輸入已知參數L，Q，Mo，E，J后，先求兩鉸鏈的支撐反力Na和Nb.
L=2; q=1000;  Mo=900;  E=200e9;  J=2e-6;
Na =(3*q*L^2/8-Mo)/L; Nb = (q*L^2/8+Mo)/L;	% 求支撐反力
x = linspace(0,L,101); dx = L/100;
M1 = Na*x(1:51)-q*x(1:51).^2/2;	% 用數組分段列出M的表達式
M2 = Nb*(L-x(52:76)) - Mo;
M3 = Nb*(L-x(77:101));
M = [M1,M2,M3];		% 列寫完整的M數組
A0 = cumtrapz(M)*dx/(E*J);	% 由M積分求轉角（未計積分常數）
Y0 = cumtrapz(A0)*dx; 	% 由轉角積分求撓度（未計積分常數）
C=[0,1;L,1]\[-Y0(1);-Y0(101)];	% 由邊界條件求積分常數Ca,Cy
Ca = C(1), Cy = C(2),
A = A0+Ca; Y=Y0+Ca*x+Cy;	% 求出轉角與撓度的完整值
subplot(3,1,1),plot(x,M),grid		%繪圖
subplot(3,1,2),plot(x,A),grid
subplot(3,1,3),plot(x,Y),grid
set(gcf,'color','w')