clear all % We need to define what is symbolic. syms a1 b1 b2 c1 c2 c3 x L EI k P real; % Our one, two and three term approximations are: w1=a1*(1-x/L)*(x/L)^2; w2=b1*(1-x/L)*(x/L)^2+b2*(1-x/L)*(x/L)^3; w3=c1*(1-x/L)*(x/L)^2+c2*(1-x/L)*(x/L)^3+c3*(1-x/L)*(x/L)^4; % This w^2, where w is uour approximation for each of the 3 cases. w1_squared=simplify(expand(w1*w1)); w2_squared=simplify(expand(w2*w2)); w3_squared=simplify(expand(w3*w3)); % These are the derivatives of w1,w2 and w3 w.r.t. x diff_w1=diff(w1,x); diff_w2=diff(w2,x); diff_w3=diff(w3,x); % These are the 2nd derivatives of w1,w2 and w3 w.r.t. x ddiff_w1=diff(diff_w1,x); ddiff_w2=diff(diff_w2,x); ddiff_w3=diff(diff_w3,x); % this is (w")^2 ddiff_w1_squared=simplify(expand(ddiff_w1*ddiff_w1)); ddiff_w2_squared=simplify(expand(ddiff_w2*ddiff_w2)); ddiff_w3_squared=simplify(expand(ddiff_w3*ddiff_w3)); % Internal energy, U, and Work, W, expresions for each case. % One term approximation. U1=EI*int(ddiff_w1_squared,x,0,L)/2+k*int(w1_squared,x,0,L)/2; W1=int(-P*w1,x,0,L); % Two term approximation. U2=EI*int(ddiff_w2_squared,x,0,L)/2+k*int(w2_squared,x,0,L)/2; W2=int(-P*w2,x,0,L); % Three term approximation. U3=EI*int(ddiff_w3_squared,x,0,L)/2+k*int(w3_squared,x,0,L)/2; W3=int(-P*w3,x,0,L); % dU/da(i) dU1_da1=diff(U1,a1); dU2_db1=diff(U2,b1); dU2_db2=diff(U2,b2); dU3_dc1=diff(U3,c1); dU3_dc2=diff(U3,c2); dU3_dc3=diff(U3,c3); % dW/da(i). These I move to right side of the eqn's b/c they are % independent of my unknowns a1, b1, b2, c1, c2, c3. B1(1,1)=diff(W1,a1); B2(1,1)=diff(W2,b1); B2(2,1)=diff(W2,b2); B3(1,1)=diff(W3,c1); B3(2,1)=diff(W3,c2); B3(3,1)=diff(W3,c3); % Now I can solve a1 (1 term approx.). a1=1; % By setting a1=1 I can seperate the coefficients of a1 from % the rest of the eqn. A1(1,1)=subs(dU1_da1); a=simplify(expand(inv(A1)*B1)); a1=a %------------------------------------------------------------------ % Result: % % a1 = % % -35/4*L^4*P/(420*EI+k*L^4) %------------------------------------------------------------------ % Now I seperate the coefficients of b1 in the two eqn's I have % (for b1, b2) and I store them in the first column of matrix A2. b1=1;b2=0; A2(1,1)=subs(dU2_db1); A2(2,1)=subs(dU2_db2); % Now I seperate the multipliers of b2 in the two eqn's I have % (for b1, b2) and I store them in the second column of matrix A2. b1=0;b2=1; A2(1,2)=subs(dU2_db1); A2(2,2)=subs(dU2_db2); % Solution = inv(A2)*B2 b=simplify(expand(inv(A2)*B2)); b1=b(1,1) b2=b(2,1) %------------------------------------------------------------------ % Result: % % b1 = % % -70*L^4*P*(6048*EI+k*L^4)/(6773760*EI^2+29568*EI*k*L^4+5*k^2*L^8) % % % b2 = % % 42*L^4*P*(6720*EI+k*L^4)/(6773760*EI^2+29568*EI*k*L^4+5*k^2*L^8) % %------------------------------------------------------------------ % Same procedure for c1, c2, c3. c1=1;c2=0;c3=0; A3(1,1)=subs(dU3_dc1); A3(2,1)=subs(dU3_dc2); A3(3,1)=subs(dU3_dc3); c1=0;c2=1;c3=0; A3(1,2)=subs(dU3_dc1); A3(2,2)=subs(dU3_dc2); A3(3,2)=subs(dU3_dc3); c1=0;c2=0;c3=1; A3(1,3)=subs(dU3_dc1); A3(2,3)=subs(dU3_dc2); A3(3,3)=subs(dU3_dc3); c=simplify(expand(inv(A3)*B3)); c1=c(1,1) c2=c(2,1) c3=c(3,1) %------------------------------------------------------------------ % Result: % % c1 = % % -105/2*L^4*P*(23950080*EI^2+ % +15552*EI*k*L^4+k^2*L^8)/(20118067200*EI^3+ % +92897280*EI^2*k*L^4+35640*EI*k^2*L^8+k^3*L^12) % % % c2 = % % 21*L^4*P*(39916800*EI^2+57600*EI*k*L^4+ % +7*k^2*L^8)/(20118067200*EI^3+ % +92897280*EI^2*k*L^4+35640*EI*k^2*L^8+k^3*L^12) % % % c3 = % % -231/2*L^8*P*k*(5040*EI+L^4*k)/(20118067200*EI^3+ % +92897280*EI^2*k*L^4+ % +35640*EI*k^2*L^8+k^3*L^12) %------------------------------------------------------------------ % I substitute a1, b1, b2, c1, c2, c3 in w1, w2, w3 respectively. w1sym=simplify(expand(subs(w1))) w2sym=simplify(expand(subs(w2))) w3sym=simplify(expand(subs(w3))) %------------------------------------------------------------------ % the approximations look like this: % % w1sym = % % 35/4*L*P*x^2*(-L+x)/(420*EI+k*L^4) % %------------------------------------------------------------------ % % w2sym = % % -14*P*x^2*(30240*EI*L^2-50400*L*x*EI+5*L^6*k-8*L^5*x*k+20160*x^2*EI+ % +3*L^4*x^2*k)/(6773760*EI^2+29568*EI*k*L^4+5*k^2*L^8) % % %------------------------------------------------------------------ % % w3sym = % % 21/2*P*x^2*(-119750400*L^2*EI^2+199584000*L*x*EI^2-77760*L^6*EI*k+ % +192960*L^5*x*EI*k-5*L^10*k^2+19*L^9*x*k^2-79833600*x^2*EI^2- % -170640*L^4*x^2*EI*k-25*L^8*x^2*k^2+55440*L^3*k*x^3*EI+ % +11*L^7*k^2*x^3)/(20118067200*EI^3+92897280*EI^2*k*L^4+ % +35640*EI*k^2*L^8+k^3*L^12) % %------------------------------------------------------------------ % this is the exact solution for the same beam without the spring % bed. (Timoshenko, Mechanics of Materials, page 685). w=-P*x^2*(3*L^2-5*L*x+2*x^2)/(48*EI); % this is the DE that gives you the exact solution for this problem % but Matlab crashes when asked to print it. % w4=dsolve('D4y=-(k/E/I)*y+P/E/I','y(0)=0,y(L)=0,Dy(0)=0','x'); % The numerical values for EI, L, k and P are: EI=1050;L=5;k=2;P=50; % We substitute them in our approximations. w=simplify(expand(subs(w))); w1=simplify(expand(subs(w1sym))); w2=simplify(expand(subs(w2sym))); w3=simplify(expand(subs(w3sym))); % The shape approximations look like this for EI=1050, L=5, k=2, P=50. %----------------------------------------------------------------------- % pretty(w1) % % 175 2 35 3 % - ---- x~ + ---- x~ % 7076 7076 %----------------------------------------------------------------------- % % pretty(w2) % % 44461550 2 14820400 3 1185618 4 % - --------- x~ + --------- x~ - --------- x~ % 600550897 600550897 600550897 % %----------------------------------------------------------------------- % pretty(w3) % % 1109865824025 2 370095011475 3 29648824149 4 % - -------------- x~ + -------------- x~ - -------------- x~ % 14987046451346 14987046451346 14987046451346 % % 4890963 5 % + -------------- x~ % 14987046451346 %----------------------------------------------------------------------- % »pretty(w) % % 25 2 25 3 4 % - --- x~ + ---- x~ - 1/504 x~ % 336 1008 %----------------------------------------------------------------------- % And now we can plot w, w1, w2, w3. figure(1) ezplot(w1,0,5) hold on ezplot(w2,0,5) hold on ezplot(w3,0,5) hold off grid axis([0 5 -.17 .01]) orient tall ylabel('Y-Direction: Deflections w1, w2, w3 [m]') xlabel('X-Direction: Distance along beam span [m]') title(' PLots of the approximated beam deflection, w1, w2, w3') % Now I set k=0 in my approximations to see if my plots for % the 2 and 3 term approximations will match the exact solution, % w,of the same problem without the the spring bed (Timoshenko, % Mechanics of Materials, page 685). EI=1050;L=5;k=0;P=50; w1=simplify(expand(subs(w1sym))); w2=simplify(expand(subs(w2sym))); w3=simplify(expand(subs(w3sym))); % Now I set k=0 and I get w2 and w3 to match the exact solution w. %----------------------------------------------------------------------- % pretty(w2) % % 25 2 25 3 4 % - --- x~ + ---- x~ - 1/504 x~ % 336 1008 %----------------------------------------------------------------------- % pretty(w3) % % 25 2 25 3 4 % - --- x~ + ---- x~ - 1/504 x~ % 336 1008 %----------------------------------------------------------------------- % pretty(w) % % 25 2 25 3 4 % - --- x~ + ---- x~ - 1/504 x~ % 336 1008 %----------------------------------------------------------------------- figure(2) ezplot(w,0,5) hold on ezplot(w2,0,5) hold on ezplot(w3,0,5) hold off grid axis([0 5 -.17 .01]) orient tall ylabel('Y-Direction: Deflections w, w2, w3 [m]') xlabel('X-Direction: Distance along beam span [m]') title(' Exact solution, Two and Three term Approximations, (w, w2 and w3), with k=0')