function cantvib3
% Last modified 7/15/2008
% This function illustrates the free vibration of a cantilever beam with
% an initial deflection curve that is both quadratic and cubic in the 
% spatial variable. This amounts to both an initial transverse load and 
% tip bending moment released at t=0.The solution is accomplished by 
% discretizing the fourth derivative of the Bernoulli-Euler beam equation
% a sixteen nodes. This script requires another function appended to the 
% end of the main script called beamrhs which calculates the righthand
% side of the differential equation for ODE45.
% This m-file was written at the University of Wyoming in
% the Electrical and Computer Engineering Department and 
% is to be distributed without cost.
clc
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
n=16;% number of nodes
M=[7 -4 1 zeros(1,n-3)];
M=[M;-4 6 -4 1 zeros(1,n-4)];
for i=3:n-2
	M=[M;zeros(1,i-3) 1 -4 6 -4 1 zeros(1,n-i-2)];
end
M=[M;zeros(1,n-4) 1 -4 5 -2];
M=[M;zeros(1,n-3) 2 -4 2];
M=((n^4)/12.3263)*M;
global A;
xoverl=[0:1/n:1];
xoverl1=xoverl(:,2:n+1);
w0=(.667*(xoverl1.^2)+.333*(xoverl1.^3));
lim1=max(w0);
x0=[w0 zeros(1,n)];
A=[zeros(n,n) eye(n,n);-M zeros(n,n)];
tspan=linspace(0,4*pi,601);
omega1t=linspace(0,4*pi,601);
[t,x]=ode45(@beamrhs,tspan,x0);
w=x(:,1:n);
sz=size(t);
w=[zeros(sz(1,1),1) w];

figure(1);clf;% Plot y(x,t) vs x/L for various times
h1=axes('YDir','reverse');
axis([-.2 1.2 -1.2*lim1 1.2*lim1])
hold on
text(0,-1,'y(x,0)=y_0[0.667(x/L)^2+0.333(x/L)^3]')
hold on
xp=[ 0  0 -.15 -.13 -.16 -.14 -.15 0];
yp=[-.37 .37 .2 .05 -.15 -.2 -.3 -.37];
patch(xp,yp,'r');
hold on
plot([-.1 0],[0 0],'k','LineWidth',[3.5])
hold on
plot([0 1.1],[0 0])
hold on
text(1.1,-.05,'x/L')
hold on
for i=1:10:150
   plot(xoverl,w(i,:),'k','LineWidth',[3.5]);
   hold on
end
xlabel('Dimensionless Distance, x/L')
ylabel('Dimensionless Displacement, y(x,t)/y_0')
box on
hold off
text(.2,1,'Press Enter to Continue');
pause


figure(2);clf;%plot y(xi,t) vs omega1t for the nodes xi
h1=axes('YDir','reverse');
axis([0 4*pi -1.2*lim1 1.2*lim1]);
hold on
for i=2:2:n
   plot(omega1t,x(:,i))
   hold on
   plot([2*pi,8],[x(301,i) .1+(i/n)],'k')
   hold on
   text(8.1,.1+(i/n),['x/L=' num2str(xoverl(1,i+1))]);
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Displacement, y(x,t)/y_0')
text(1,1,'Press Enter to Continue');
%box on
hold off
pause
[T,X]=meshgrid(omega1t,xoverl);

figure(3);clf;
axis([0 4*pi 0 1 -1.2*lim1 1.2*lim1]);
hold on
mesh(T,X,w')
colormap([0 0 0]);
view(60,30)
xlabel('Dimensionless Time, \omega_1t','rotation',-31)
ylabel('Dimensionless Location, x/L','rotation',12)
zlabel('Displacement, y(x,t)/y_0')
grid on
hold off
pause

figure(4);clf;%Animate y(x,omega1t)
h1=axes('YDir','reverse');
axis([-.2 1.2 -1.2*lim1 1.2*lim1])
hold on
xp=[ 0  0 -.15 -.13 -.16 -.14 -.15 0];
yp=[-.37 .37 .2 .05 -.15 -.2 -.3 -.37];
patch(xp,yp,'r');
hold on
plot([-.1 0],[0 0],'k','LineWidth',[3.5])
hold on
plot([0 1.1],[0 0])
L=plot(xoverl,zeros(1,n+1),'k','EraseMode','xor','LineWidth',[3.5]);
hold on
text(1.1,-.05,'x/L')
hold on
box on
text(0,-1,'y(x,0)=y_0[0.667(x/L)^2+0.333(x/L)^3]')
thandl=text(0.2,1,'Press Enter to Set Initial Condition');
xlabel('Dimensionless Distance, x/L')
ylabel('Dimensionless Displacement, y(x,t)/y_0')
pause
for i=1:sz(1,1)
   set(L,'YData',w(i,:));
   if i==1
      set(thandl,'String','Press Enter to Animate'); 
      pause
   end
   drawnow;
   pause(.01)
end
hold off
set(thandl,'String','Press Enter to Continue'); 
pause


function xdot=beamrhs(t,x)
global A;
xdot=A*x;