clear

%  求解單元 1 的剛度矩陣
x=sym ('x');        L1=sym ('L1');      % 定義x和L1為符號變量
E=3e11;        I=5.2e-7;
N=[ 1 - 3 * (x^2) / (L1^2) + 2 * (x^3) / (L1^3),  x - 2 * (x^2) / L1 + (x^3) / (L1^2),....
       3 * (x^2) / (L1^2) - 2 * (x^3) / (L1^3),  -(x^2) / L1 + (x^3) / (L1^2) ];    % 形函數
Ni=diff ( N, x, 2 ); % 對形函數進行2次微分
Nt=transpose ( Ni);
kk1=Nt * Ni;
K1=E*I*int(kk1,0,'L1');      % 對kk1進行積分，積分區間0－L1
K1=subs (K1, 'L1', 5);         % 將K1矩陣中的L1代換成5，重新計算K1


%  求解單元 2 的剛度矩陣
x=sym ('x');        L2=sym ('L2');    % 定義x和L2為符號變量
E=3e11;        I=5.2e-7;
N=[ 1 - 3 * (x^2) / (L2^2) + 2 * (x^3) / (L2^3),  x - 2 * (x^2) / L2 + (x^3) / (L2^2),....
      3 * (x^2) / (L2^2) - 2 * (x^3) / (L2^3),  -(x^2) / L2 + (x^3) / (L2^2) ];     % 形函數
Ni=diff ( N, x, 2 );     % 對形函數進行2次微分
Nt=transpose ( Ni);
kk2=Nt * Ni;
K2=E* I * int (kk2, 0, 'L2');     % 對kk2進行積分，積分區間0－L2
K2=subs (K2, 'L2', 5);            % 將K1矩陣中的L1代換成5，重新計算K1


%  合并剛度矩陣
K1=[K1,zeros(4,2);zeros(2,6)];
K2=[zeros(2,6);zeros(4,2),K2];
KK=K1 + K2;


%  求解最后結果
m=KK(3:6,3:6);        % 將KK矩陣的3-6行和3-6列取出，組成一個4X4的矩陣
format  rat           % 將計算結果以分數形式表示
mm= inv ( m );
F=[0;  0;  -100;  0];
format  long          % 將計算結果以16位有效數字表示
u=mm * F