Matlab Beam and Differential Equation Tools

Edit: This post and the code have been getting quite a few hits.  If you download and use/adapt these Matlab scripts, could you comment to let me know how they are being used?

Beam

I have been playing around with using Matlab to make learning engineering concepts more interactive for a few years now and have developed a couple of programs that are both fun and educational.  One is designed to make learning how to draw beam bending moment diagrams easier and the other to visualise what the solutions to differential equations tell us. I have made them freely available here

Bending Moment Diagram Tool

This tool is designed to generate examples of bending moment (and shear force) diagram problems.  Each problem is generated randomly so it is possible to practice on many sample problems.  A choice of what type of load (point, UDL, triangle, point moment) is possible and it is also possible to change the sign conventions from the defaults, just to make everything more confusing!.  The code is below – it should run from Matlab if expanded and saved to a .m file

%MAIN FUNCTION WHICH LAYS OUT WINDOW AND RESPONDS TO INPUT
function beam4()
 
% This function specifies the GUI for producing a bending moment diagram
% problem.
PointYN=0;
UDLYN=0;
TriangYN=0;
PointMomYN=0;
length=0;supp1=0 ; supp2=0; NumPload=0; NumDload=0; NumTload=0; NumPin=0; NumPMoment=0;...
 Pload=0; UDL=0; Tload=0; Pin=0; Mom=0;Totload=0;Reaction(2)=0;LoadMat=0; Shear=0;BM=0;
 
%Set default sign conventions
BMSign=-1; %sagging positive
ShearSign=-1; %clockwise positive
 
% SPECIFY ENTIRE WINDOW FOR GUI
fh=figure('name', 'Bending Moment Diagram Tutorial','MenuBar',...
 'none','Units','normalized','color',[0.95 0.95 0.95],...
 'Position', [0.05,0.05,0.9,0.7]);
 
%SPECIFY PROBLEM DEFINITION BOX AND DRAWING OF PROBLEM
%Text for check boxes that select types of load
inpanel = uipanel('Parent',fh,'Title','Problem to include',...
 'Position',[0.02 0.72 0.22 0.22],'fontsize',12,...
 'FontWeight','bold');
Points=uicontrol(inpanel,'Style','text',...
 'String','Point loads','Units','normalized',...
 'Position',[0,0.8,0.4,0.1],'HorizontalAlignment','left');
UDLs=uicontrol(inpanel,'Style','text',...
 'String','Uniform loads','Units','normalized', 'Position',[0,0.6,0.4,0.1],'HorizontalAlignment','left');
Triangs=uicontrol(inpanel,'Style','text',...
 'String','Triangular loads','Units','normalized','Position',[0,0.4,0.4,0.1],'HorizontalAlignment','left');
PointMoms=uicontrol(inpanel,'Style','text',...
 'String','Point moments','Units','normalized','Position',[0,0.2,0.4,0.1],'HorizontalAlignment','left');
 
%Check boxes for types of load
PointsCh = uicontrol(inpanel,'Style','CheckBox',...
 'Units','normalized','Callback',@getPointsCh,...
 'Position',[0.4,0.8,0.4,0.1],'HorizontalAlignment','left');
UDLSCh = uicontrol(inpanel,'Style','CheckBox',...
 'Units','normalized','Callback',@getUDLsCh,...
 'Position',[0.4,0.6,0.4,0.1],'HorizontalAlignment','left');
TraingsCh = uicontrol(inpanel,'Style','CheckBox',...
 'Units','normalized','Callback',@getTriangsCh,...
 'Position',[0.4,0.4,0.4,0.1],'HorizontalAlignment','left');
PointMomsCh = uicontrol(inpanel,'Style','CheckBox',...
 'Units','normalized','Callback',@getPointMomsCh,...
 'Position',[0.4,0.2,0.4,0.1],'HorizontalAlignment','left');
 
%Call back functions for check boxes
 function getPointsCh(hObject,eventdata)
 PointYN=get(hObject,'value');
 end
 
function getUDLsCh(hObject,eventdata)
 UDLYN=get(hObject,'value');
 end
 
function getTriangsCh(hObject,eventdata)
 TriangYN=get(hObject,'value');
 end
 
function getPointMomsCh(hObject,eventdata)
 PointMomYN=get(hObject,'value');
 end
 
%SPECIFY REACTION FORCE BOX AND SFD PLOTTING
inpanel2 = uipanel('Parent',fh,'Title','Reactions and SFD',...
 'Position',[0.02 0.42 0.22 0.22],'fontsize',12,...
 'FontWeight','bold');
 
%SPECIFY BMD BOX
inpanel3 = uipanel('Parent',fh,'Title','BMD',...
 'Position',[0.02 0.12 0.22 0.22],'fontsize',12,...
 'FontWeight','bold');
 
%SPECIFY SIGN CONVENTION BOX
inpanel4 = uipanel('Parent',fh,'Title','Sign Conventions',...
 'Position',[0.02 0.01 0.22 0.09],'fontsize',8,...
'FontWeight','bold');
 
%Text for bending moments
uicontrol(inpanel4,'Style','text',...
 'String','BM','fontweight','bold','Units','normalized','fontsize',8,...
 'Position',[0,0,0.4,0.4],'HorizontalAlignment','left');
uicontrol(inpanel4,'Style','text',...
 'String','Sagging pos','Units','normalized','fontsize',8,...
 'Position',[0.15,0,0.3,0.4],'HorizontalAlignment','left');
uicontrol(inpanel4,'Style','text',...
 'String','Sagging neg','Units','normalized','fontsize',8,...
 'Position',[0.59,0,0.4,0.4],'HorizontalAlignment','left');
 
%Text for shear force
uicontrol(inpanel4,'Style','text',...
 'String','Shear','fontweight','bold','Units','normalized','fontsize',8,...
 'Position',[0,0.45,0.4,0.4],'HorizontalAlignment','left');
uicontrol(inpanel4,'Style','text',...
 'String','Clockwise pos','Units','normalized','fontsize',8,...
 'Position',[0.15,0.45,0.3,0.4],'HorizontalAlignment','left');
uicontrol(inpanel4,'Style','text',...
 'String','Clockwise neg','Units','normalized','fontsize',8,...
 'Position',[0.59,0.45,0.4,0.4],'HorizontalAlignment','left');
%Checkbox callbacks
CHB1=uicontrol(inpanel4,'Style','CheckBox','Value',1,...
 'Units','normalized','Callback',@getSign1Ch,...
 'Position',[0.45,0.12,0.1,0.2],'HorizontalAlignment','left');
CHB2=uicontrol(inpanel4,'Style','CheckBox',...
 'Units','normalized','Callback',@getSign2Ch,...
 'Position',[0.89,0.12,0.1,0.2],'HorizontalAlignment','left');
 
CHB3=uicontrol(inpanel4,'Style','CheckBox','Value',1,...
 'Units','normalized','Callback',@getSign3Ch,...
 'Position',[0.45,0.59,0.1,0.2],'HorizontalAlignment','left');
 
CHB4=uicontrol(inpanel4,'Style','CheckBox',...
 'Units','normalized','Callback',@getSign4Ch,...
 'Position',[0.89,0.59,0.1,0.2],'HorizontalAlignment','left');
 
%Callback functions that setup sign conventions
function getSign1Ch(hObject,eventdata)
 BMSign=get(hObject,'value');
 if BMSign==1
 set(CHB2,'Value',0)
 BMSign=-1;
% else
% set(CHB2,'Value',1)
% BMSign=-1;
 end
end
function getSign2Ch(hObject,eventdata)
 BMSign=get(hObject,'value');
 if BMSign==1
 set(CHB1,'Value',0)
 BMSign=1;
% else
% set(CHB1,'Value',1)
% BMSign=1;
 end
end
 
function getSign3Ch(hObject,eventdata)
 ShearSign=get(hObject,'value');
 if ShearSign==1
 set(CHB4,'Value',0)
 ShearSign=-1;
% else
% set(CHB4,'Value',1)
% ShearSign=1;
 end
end
 
function getSign4Ch(hObject,eventdata)
 ShearSign=get(hObject,'value');
 if ShearSign==1
 set(CHB3,'Value',0)
 ShearSign=1;
% else
% set(CHB3,'Value',1)
% ShearSign=-1;
 end
end
 
%BUTTON TO DRAW PROBLEM
uicontrol(inpanel,'Style','pushbutton','String','Draw Problem',...
 'Callback', {@stage1}...
 ,'Units','normalized','Position',[0.55,0.4,0.4,0.25]);
 
%BUTTON TO CALC REACTIONS
uicontrol(inpanel2,'Style','pushbutton','String','Show Reactions',...
 'Callback', {@stage2},...
 'Units','normalized','Position',[0.55,0.6,0.4,0.25]);
 
%BUTTON TO CALC SFD
uicontrol(inpanel2,'Style','pushbutton','String','Draw SFD',...
 'Callback', @stage3,...
 'Units','normalized','Position',[0.55,0.2,0.4,0.25]);
%BUTTON TO CALC BMD
uicontrol(inpanel3,'Style','pushbutton','String','Draw BMD',...
 'Callback', @stage4,...
 'Units','normalized','Position',[0.55,0.4,0.4,0.25]);
%functions that happens when define problem button is pressed
 function []=stage1(hObject,eventdata)
 
 %call params function which sets up problem with random parameters
 %for each type of load selected in the gui
 [length supp1 supp2 NumPload NumDload NumTload NumPin NumPMoment Pload UDL Tload...
 Pin Mom]=params(PointYN,UDLYN,TriangYN,PointMomYN);
 
 %plots problem using the above generated parameters
 plotting1(length, supp1, supp2, NumPload, NumDload, NumTload, NumPin, NumPMoment,...
 Pload, UDL, Tload, Pin, Mom)
 
 %puts all loads in one UDL array
 [LoadMat]=LM(NumPload, NumDload, NumTload, ...
 Pload, UDL, Tload,length);
 end
 
%functions than happens when the show reactions button is pressedc
 function stage2(hobject,eventdata)
 
 [Totload, Reaction Shear BM]=react(length,supp1,supp2, LoadMat, NumPMoment, Mom,BMSign, ShearSign);
 
 plotting2(Shear, Reaction,length, supp1, supp2,NumPin,Pin)
 end
 
%function that happens when draw SFD button is pressed
 function stage3(hobject,eventdata)
 plotting3(Shear)
 end
 
%function that happens when draw BM button is pressed
 function stage4(hobject,eventdata)
 plotting4(BM, length)
 end
 
%Development details
devpanel = uipanel('Parent',fh,'Position',[0.5 0.01 0.3 0.05]);
t(1)={'Developed by and copyright of Dr Martin Gillie, University of Manchester'};
t(2)={'Comments to martin.gillie@manchester.ac.uk'};
uicontrol(devpanel,'Style','text',...
 'String',t,'Units','normalized',...
 'HorizontalAlignment','center','Position',[0.02 0.02 0.9 0.9]);
end
 
%FUNCTION SETS RANDOM PARAMETERS FOR PROBLEM
function [length supp1 supp2 NumPload NumDload NumTload NumPin NumPMoment Pload UDL...
 Tload Pin Mom]=params(PointYN,UDLYN,TriangYN,PointMomYN)
 
NumPload=0;NumDload=0;NumTload=0;NumPin=0;NumPMoment=0;Pload=0;
UDL=0;Tload=0;Pin=0;Mom=0;
 
%set beam parameters
maxlength=20; %max allowed beam length
minlength=6; %min allowed beam length
length=random('unif',minlength,maxlength); %length of beam
supp1=(random('unif',0,length/3)); %position of 1st support
supp2=(random('unif',length/3+2,length)); %position of 2nd support
 
for i=0:maxlength
 if length >=i && length < i+1
 length=i; %make length a round number but not integer
 end
 if supp1>=i && supp1<i+1
 supp1=i; %make length a round number but not integer
 end
 if supp2>=i && supp2<i+1
 supp2=i; %make length a round number but not integer
 end
end
 
if PointYN==1
 NumPload=int8(random('unif',1,3)); %number of point loads
 for i=1:NumPload;
 Pload(i,1)=cast(int8(random('unif',-50,50)),'single'); %magnitude
 if Pload(i,1)==0
 Pload(i,1)=1; %ensure no zero loads
 end
 Pload(i,2)=cast(int8(random('unif',0,length)),'single'); %position
 if Pload(i,2)==supp1
 Pload(i,2)=Pload(i,2)+1; %no loads at supports
 end
 if Pload(i,2)==supp2
 Pload(i,2)=Pload(i,2)+1; %no loads at supports
 end
 end
end
 
if UDLYN==1
NumDload=int8(random('unif',1,2)); %number of UDLs
 for i=1:NumDload
 if NumDload==1
 UDL(i,1)=cast(int8(random('unif',0,length-2)),'single'); %any length for single load
 UDL(i,2)=cast(int8(random('unif',UDL(i,1)+1,length)),'single');
 else
 if i==1
 UDL(i,1)=cast(int8(random('unif',0,length/4)),'single'); %max start at l/4 for 1st load
 UDL(i,2)=cast(int8(random('unif',UDL(i,1)+1,length/2)),'single');
 else
 UDL(i,1)=cast(int8(random('unif',length/2,length-1)),'single'); %min start at l/2 for 2nd load
 UDL(i,2)=cast(int8(random('unif',UDL(i,1)+1,length)),'single');
 end
 end
 UDL(i,3)=cast(int8(random('unif',-10,10)),'single');
 if UDL(i,3)==0
 UDL(i,3)=1; %ensure no zero loads
 end
 
 end
end
 
if TriangYN==1
 NumTload=int8(random('unif',1,1)); %number of triang loads
 if NumTload>0 %assign properties to tloads
 for i=1:NumTload
 Tload(1)=cast(int8(random('unif',0, length-2)),'single'); %start
 Tload(2)=cast(int8(random('unif',Tload(1)+1, length)),'single'); %stop
 Tload(3)=cast(int8(random('unif',1,5)),'single'); %magnitude
 end
 end
end
 
if PointMomYN==1
 NumPMoment=int8(random('unif',1,2)); %number of point moments
 if NumPMoment>0
 for i=1:NumPMoment
 Mom(i,1)=cast(int8(random('unif',0,length)),'single'); %location
 Mom(i,2)=cast(int8(random('unif',-200,200)),'single');
 if Mom(i,2)==0
 Mom(i,2)=1; %no zero loads
 end
 end
 end
end
end
 
%FUNCTION PLOTS PROBLEM
function []=plotting1(length, supp1, supp2, NumPload, NumDload, NumTload, NumPin,...
 NumPMoment, Pload, UDL, Tload, Pin, Mom)
 
%set up axes
subplot (3,3,[2 3])
cla
axis([0-length*0.1,length*1.1,-1,1])
set(gca,'YTick',[])
set(gca,'YTickLabel',{})
set(get(gca,'XLabel'),'String','Position (m)')
 
%draw beam
line([0 length],[0 0],'linewidth',5);
%draw point loads, if any
if NumPload>0
 for i=1:NumPload;
 APos(1)=0.5; %coord for drawing arrows
 APos(2)=0.15; %coord for drawinfg arrow heads
 APos(3)=0.045;
 APos(4)=0.6; %text position
 if Pload(i,1)<0
 APos=-APos; %draw negative loads upwards
 end
 
line([Pload(i,2) Pload(i,2)],[APos(1),APos(3)],'linewidth',2,'color','red')
 line([Pload(i,2)-length/70 Pload(i,2)],[APos(2),APos(3)],'linewidth',2,'color','red')
 line([Pload(i,2)+length/70 Pload(i,2)],[APos(2),APos(3)],'linewidth',2,'color','red')
 
 text('Position',[Pload(i,2)-length/50 APos(4)],'String',[num2str(Pload(i,1),3),'kN'])
 end
end
 
%draw UDLs, if any
if NumDload>0
 for i=1:NumDload;
 if UDL(i,3)>0 %plot positve UDLs above axis
 UDL(i,1)
 rectangle('Position',[UDL(i,1) 0.03 UDL(i,2)-UDL(i,1) 0.09],'FaceColor','red')
 text('Position',[UDL(i,1) 0.27],'String',[num2str(UDL(i,3),3),'kN/m'])
 else %and negative ones below
 rectangle('Position',[UDL(i,1) -0.12 UDL(i,2)-UDL(i,1) 0.09],'FaceColor','red')
 text('Position',[UDL(i,1) -0.27],'String',[num2str(UDL(i,3),3),'kN/m'])
 end
 end
end
 
%draw Tloads, if any
if NumTload>0
hold all
 for i=1:NumTload
 plot([Tload(1) Tload(2) Tload(2)],[0.045 0.8 0.045],'linewidth',2,'color','red')
 text('Position',[Tload(2)+length/80 0.8],'String',[num2str(Tload(3),3) 'kN peak'])
 end
end
%draw pins, if any
 if NumPin>0
 for i=1:NumPin
 hold all
 plot(Pin(1), 0,'Marker','o','MarkerFaceColor','w','MarkerEdgeColor','k')
 end
 end
 
 %draw point moments if any
 if NumPMoment>0
 for i=1:NumPMoment
 hold all
 plot(Mom(i,1), 0,'Marker','o','MarkerFaceColor','r','MarkerEdgeColor','r')
 text('Position',[Mom(i,1)-length/20 -0.15],'String',[num2str(Mom(i,2),3), 'kNm'])
 end
 end
 
%draw supports
 patch([supp1-length/60 supp1 supp1+length/60],[-0.2 -0.03 -0.2],[0.1 0.1 0.1]);
 patch([supp2-length/60 supp2 supp2+length/60],[-0.2 -0.03 -0.2],[0.1 0.1 0.1]);
 
end
 
%FUNCTION CONVERTS ALL LOADS TO UDLS AND PLACES IN AN ARRAY
function[LoadMat]=LM(NumPload, NumDload, NumTload, ...
 Pload, UDL, Tload,length)
 
LoadMat=zeros(NumPload+NumDload*99+NumTload*999,3);
Row=0; %number of UDLs in total
if NumPload >0
 for i=1:NumPload;
 Row=Row+1;
 LoadMat(Row,1)=Pload(i,2); %set start position of UDL
 LoadMat(Row,2)=Pload(i,2)+length/1000; %set end position of UDL
 LoadMat(Row,3)=Pload(i,1)/(LoadMat(i,2)-LoadMat(i,1));
 end
end
 
if NumDload >0
 for i=1:NumDload;
 UDLSize=(UDL(i,2)-UDL(i,1))/100;
 for j=0:99
 Row=Row+1;
 LoadMat(Row,1)=UDL(i,1)+UDLSize*j; %set start position of UDL
 LoadMat(Row,2)=LoadMat(Row,1)+UDLSize; %set end position of UDL
 LoadMat(Row,3)=UDL(i,3); %set magnitude of UDL
 end
 end
end
 
if NumTload >0
 for i=1:NumTload;
 UDLSize=(Tload(i,2)-Tload(i,1))/1000;
 for j=0:999
 Row=Row+1;
 LoadMat(Row,1)=Tload(i,1)+UDLSize*j; %set start position of UDL
 LoadMat(Row,2)=LoadMat(Row,1)+UDLSize; %set end position of UDL
 LoadMat(Row,3)=(Tload(i,3)/(Tload(i,2)-Tload(i,1)))*(LoadMat(Row,1)-Tload(i,1)); %set magnitude of UDL
 end
 end
end
end
 
%FUNCTION CALCULATES TOTAL LOAD, SUPPORT REACTIONS, SF DIST AND BM DIST
function [Totload Reaction Shear BM]=react(length,supp1,supp2, LoadMat, NumPMom, Mom,BMSign,ShearSign)
 
Reaction=0;Totload=0;MomSupp1=0;IncMom=0;ForcePos=0;Shear=0;Mom;
rows=size(LoadMat,1); % Used in calcs
 
%CALCULATE TOTAL LOAD AND SUPPORT REACTIONS
%obtain total load, and total moment about support 1
if size(LoadMat,2) ~=1 %Do not enter if no loads
 for i=1:size(LoadMat,1)
 Totload=Totload+(LoadMat(i,3)*(LoadMat(i,2)-LoadMat(i,1))); %Total load
 ForcePos=(LoadMat(i,2)+LoadMat(i,1))/2; %Midpoint of each UDL load
 IncMom=LoadMat(i,3)*(LoadMat(i,2)-LoadMat(i,1))*(ForcePos-supp1); %Moment due to each UDL about support 1
 MomSupp1=MomSupp1+IncMom; %Total Mom around support 1
 end
end
 
%account for effect of all point moments on overall moment of beam
for i=1:NumPMom
 MomSupp1=Mom(i,2)+MomSupp1; %moment about supp1 due to loads and moments
end
 
%calculate reaction forces
Reaction(2)=MomSupp1/(supp2-supp1); %support reaction 2
Reaction(1)=Totload-Reaction(2); %supprot reaction 1
%CALCULATE SHEAR FORCE DISTRIBUTION
 
%Add in support reactions to load matrix
rows=rows+1;
LoadMat(rows,1)=supp1; %add support reaction 1
LoadMat(rows,2)=supp1+length/1000;
LoadMat(rows,3)=-1*Reaction(1)/(LoadMat(rows,2)-LoadMat(rows,1)); %different sign for reactions so *-1
rows=rows+1;
LoadMat(rows,1)=supp2; %add support reaction 2
LoadMat(rows,2)=supp2+length/1000;
LoadMat(rows,3)=-1*Reaction(2)/(LoadMat(rows,2)-LoadMat(rows,1)); %different sign for reactions so *-1
 
%sort load matrix to get UDL start points in increasing order.
LoadMat=sortrows(LoadMat,1);
 
%fill matrix with shear forces
Shear=zeros(1000,2);
for i=1:1000 %divide beam into 1000 chunks
 Xpos=(i-1)*length/1000;
 Shear(i,1)=Xpos; %position at which shear is being calculate
 for j=1:rows %calc shear at beginning of each chunk
 if LoadMat(j,1) < Xpos %only consider loads to left of current position
 Shear(i,2)=Shear(i,2)+((LoadMat(j,2)-LoadMat(j,1))*LoadMat(j,3));
 end
 end
 Shear(i,2)=Shear(i,2)*ShearSign; %account for shear sign convention
end
 
%CALCULATE BENDING MOMENT DISTRIBUTION
BM=zeros(1000,2);
for i=1:1000 %divide beam into 1000 chunks
 Xpos=(i-1)*length/1000;
 for j=1:rows %calc moment at beginning of each chunk
 BM(i,1)=Xpos; %position at which moment is being calculated
 if LoadMat(j,1) < Xpos %only consider loads to left of current position
 IncBM=((LoadMat(j,2)-LoadMat(j,1))*LoadMat(j,3))*(Xpos-LoadMat(j,1));
 BM(i,2)=BM(i,2)+IncBM;
 end
 end
 for k=1:NumPMom %account for point moments
 if Mom(k,1)<Xpos
 BM(i,2)=BM(i,2)-Mom(k,2);
 end
 end
 BM(i,2)=BM(i,2)*BMSign; %account for BM Sign convention
end
end
 
%FUNCTION PLOTS REACTION FORCES IN SECOND FIGURE
function []=plotting2(Shear,Reaction, length, supp1, supp2,NumPin,Pin)
 
%set up plot area
subplot (3,3,[5 6])
cla
set(gca,'YTick',[])
set(get(gca,'YLabel'),'String','')
set(get(gca,'XLabel'),'String','Position (m)')
 
%set axes limits etc
minax=min(Shear,[],1)*1.1;
maxax=max(Shear,[],1)*1.1;
 
if minax(2)==0
 minax(2)=-10;
end
if maxax(2)==0
 maxax(2)=10;
end
 
axis([0-length*0.1,length*1.1,minax(2)*1.1,maxax(2)*1.1])
 
%draw beam
line([0 length],[0 0],'linewidth',1);
 
%draw pins, if any
if NumPin>0
 for i=1:NumPin
 hold all
 plot(Pin(1), 0,'Marker','o','MarkerFaceColor','w','MarkerEdgeColor','k')
 end
end
 
%draw supports
posrange=(minax(2)*-1+maxax(2))/2;
patch([supp1-length/60 supp1 supp1+length/60],[-0.2*posrange -0.03*posrange -0.2*posrange],[0.1 0.1 0.1]);
patch([supp2-length/60 supp2 supp2+length/60],[-0.2*posrange -0.03*posrange -0.2*posrange],[0.1 0.1 0.1]);
 
text('Position', [supp1 -0.3*posrange],'String', [num2str(Reaction(1),3), 'kN']);
text('Position', [supp2 -0.3*posrange],'String', [num2str(Reaction(2),3), 'kN']);
end
 
%FUNCTION PLOTS SFD IN SECOND FIGURE
function[]=plotting3(Shear)
 
subplot (3,3,[5 6])
set(gca,'YTickMode','auto')
set(get(gca,'YLabel'),'String','')
set(get(gca,'YLabel'),'String','Shear Force (kN)')
chunks=size(Shear,1);
for i=1:chunks-1
 line([Shear(i) Shear(i+1)],[Shear(i,2) Shear(i+1,2)],'Color',[0,0,0],'linewidth',2)
end
line([Shear(chunks,1) Shear(chunks,1)],[Shear(chunks,2) 0],'Color',[0,0,0],'linewidth',2)
 
end
 
%FUNCTION PLOTS BMD IN THIRD FIGURE
function[]=plotting4(BM,length)
 
subplot (3,3,[8 9])
cla
set(get(gca,'YLabel'),'String','Bending Moment (kNm)')
set(get(gca,'XLabel'),'String','Position (m)')
%set axes limits
minax=min(BM,[],1);
maxax=max(BM,[],1);
if minax(2)>-2
 minax(2)=-2;
end
if maxax(2)<2
 maxax(2)=2;
end
 
axis([0-length*0.1,length*1.1,minax(2)*1.1,maxax(2)*1.1])
%draw beam
line([0 length],[0 0],'linewidth',1);
chunks=size(BM,1);
for i=1:chunks-1
 line([BM(i) BM(i+1)],[BM(i,2) BM(i+1,2)],'Color',[0,0,0],'linewidth',2)
end
line([BM(chunks,1) BM(chunks,1)],[BM(chunks,2) 0],'Color',[0,0,0],'linewidth',2)
 
end

Differential Equation Solver

This tool solves linear or non-linear first or second order differential equations and plots the results graphically.  It is possible to see and compare the effects of damping and various forms of driving force.  The Matlab script (which should just run from Matlab if expanded and copied to a *.m file) is below and an instruction sheet here


function diffeqnview7()
% This functions specifies the GUI for solving second order differential
%equations associated with oscillating masses etc. It is concerned only
%with the graphical interface and obtaining user input data. It calls
%function &quot;solution&quot;, below, which actually produces the mathematical
%solution. Within diffenqview2 are a number of nested functions that
%are called when the user does something, like enter a new value for a
%parameter.

% SPECIFY ENTIRE WINDOW FOR GUI
fh=figure('name', 'Single Degree of Freedom Oscillator','MenuBar',...
 'none','Units','normalized','color',[0.95 0.95 0.95],...
 'Position', [0.05,0.05,0.9,0.7]);

% CONSTRUCT INPUT PARAMETERS BIT OF GUI
%Initialize parameters to boring values for system 1
mass(1)=1;
stiffness{1}='1'; %entered as string for flexibility in equations
damping{1}='0.1'; %entered as string for flexibility in equations
time(1)=20; %Total time over which solution is calculated
ivelocity(1)=1; %Initial velocity of mass
iposition(1)=1; %Initial position of mass
dftype(1)=1; %Type of driving force - number a &quot;code&quot; for none.
dfperiod(1)=0; %Period of driving force function
dfmag(1)=0; %Peak magnitude of driving force
dfsttime(1)=0; %Start time of driving force
dfsptime(1)=0; %Stop time of driving force

%Initialize parameters to boring values for system 2. All values in
%vector format.
mass(2)=1;
stiffness{2}='1';
damping{2}='0.2';
time(2)=20; %Total time over which solution is calculated
ivelocity(2)=1; %Initial velocity of mass
iposition(2)=1; %Initial position of mass
dftype(2)=1; %Type of driving force - number a &quot;code&quot; for none.
dfperiod(2)=0; %Period of driving force function
dfmag(2)=0; %Peak magnitude of driving force
dfsttime(2)=0; %Start time of driving force
dfsptime(2)=0; %Stop time of driving force

%Input parameters panel
inpanel = uipanel('Parent',fh,'Title','Input Parameters',...
 'Position',[0.02 0.25 0.22 0.7],'fontsize',12,...
 'FontWeight','bold');

%Button to make stuff happen - solve, plot etc.
%@solve tells Matlab to call &quot;solve&quot;, a callback function, when the
%button is pressed that in turn calls the the solution function.
uicontrol(inpanel,'Style','pushbutton','String','Solve',...
 'Callback', @solve,...
 'Units','normalized','Position',[0.2,0.01,0.6,0.05]);

%Parameter name boxes. These are just identifying bits of text
%next to the boxes where user input is made

%Headings
sys1=uicontrol(inpanel,'Style','text',...
 'String','System 1','Units','normalized','fontsize',10,...
 'HorizontalAlignment','center','foregroundcolor','k',...
 'FontWeight','bold','Position',[0. 0.87 0.5 0.1]);
sys2=uicontrol(inpanel,'Style','text',...
 'String','System 2','Units','normalized','fontsize',10,...
 'HorizontalAlignment','center','foregroundcolor','r',...
 'FontWeight','bold','Position',[0.5 0.87 0.5 0.1]);
align([sys1 sys2],'Distribute','None');

%Basic input parameters labels for unforced oscillations in sys 1
massl=uicontrol(inpanel,'Style','text',...
 'String','Mass','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.85 0.4 0.05]);
stiffnessl=uicontrol(inpanel,'Style','text',...
 'String','Damping','Units','normalized','HorizontalAlignment',...
 'left','Position',[0.03 0.84 0.4 0.05]);
dampingl=uicontrol(inpanel,'Style','text',...
 'String','Stiffness','Units','normalized','HorizontalAlignment',...
 'left','Position',[0.03 0.83 0.4 0.05]);
ipositionl=uicontrol(inpanel,'Style','text',...
 'String','Initial position','Units','normalized','HorizontalAlignment',...
 'left', 'Position',[0.03 0.82 0.4 0.05]);
ivelocityl=uicontrol(inpanel,'Style','text',...
 'String','Initial velocity','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.81 0.4 0.05]);
ttimel=uicontrol(inpanel,'Style','text',...
 'String','Time period','Units','normalized',...
 'HorizontalAlignment','left', 'Position',[0.03 0.5 0.4 0.05]);
align([massl stiffnessl dampingl ipositionl ivelocityl ...
 ttimel],'None','Distribute');

 %and in sys 2
massl2=uicontrol(inpanel,'Style','text',...
 'String','Mass','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.85 0.4 0.05]);
stiffnessl2=uicontrol(inpanel,'Style','text',...
 'String','Damping','Units','normalized','HorizontalAlignment',...
 'left','Position',[0.53 0.84 0.4 0.05]);
dampingl2=uicontrol(inpanel,'Style','text',...
 'String','Stiffness','Units','normalized','HorizontalAlignment',...
 'left','Position',[0.53 0.83 0.4 0.05]);
ipositionl2=uicontrol(inpanel,'Style','text',...
 'String','Initial position','Units','normalized','HorizontalAlignment',...
 'left', 'Position',[0.53 0.82 0.4 0.05]);
ivelocityl2=uicontrol(inpanel,'Style','text',...
 'String','Initial velocity','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.81 0.4 0.05]);
ttimel2=uicontrol(inpanel,'Style','text',...
 'String','Time period','Units','normalized',...
 'HorizontalAlignment','left', 'Position',[0.53 0.5 0.4 0.05]);
align([massl2 stiffnessl2 dampingl2 ipositionl2 ivelocityl2 ...
 ttimel2],'None','Distribute');

%Additional labels for driving force terms sys1
dftypel=uicontrol(inpanel,'Style','text',...
 'String','Driving force','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.4 0.15 0.06]);
dfmagl=uicontrol(inpanel,'Style','text',...
 'String','Peak magnitude','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.2 0.2 0.06]);
dfperiodl=uicontrol(inpanel,'Style','text',...
 'String','Period (sin only)','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.25 0.2 0.06]);
dfsttimel=uicontrol(inpanel,'Style','text',...
 'String','Start time','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.15 0.2 0.05]);
dfsptimel=uicontrol(inpanel,'Style','text',...
 'String','Stop time','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.03 0.1 0.2 0.05]);
align([dftypel dfmagl dfperiodl dfsttimel dfsptimel ],'None','Distribute');

%andsys2
dftypel2=uicontrol(inpanel,'Style','text',...
 'String','Driving force','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.4 0.15 0.06]);
dfmagl2=uicontrol(inpanel,'Style','text',...
 'String','Peak magnitude','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.2 0.2 0.06]);
dfperiodl2=uicontrol(inpanel,'Style','text',...
 'String','Period (sin only)','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.25 0.2 0.06]);
dfsttimel2=uicontrol(inpanel,'Style','text',...
 'String','Start time','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.15 0.2 0.05]);
dfsptimel2=uicontrol(inpanel,'Style','text',...
 'String','Stop time','Units','normalized',...
 'HorizontalAlignment','left','Position',[0.53 0.1 0.2 0.05]);
align([dftypel2 dfmagl2 dfperiodl2 dfsttimel2 dfsptimel2],'None','Distribute');

% %Parameter value boxes - all within parameters panel
% %@get... tells Matlab what to do when data is entered
%
%Basic input parameters for unforced oscillations sys1
massb1=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback', @getmass,...
 'Position',[0.3 0.85 0.15 0.05]);
dampingb1=uicontrol(inpanel,'Style','edit',...
 'String','0.1','Units','normalized',...
 'BackgroundColor','White','Callback', @getdamping,...
 'Position',[0.3 0.57 0.15 0.05]);
stiffnessb1=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback',@getstiffness,...
 'Position',[0.3 0.56 0.15 0.05]);
ipositionb1=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback', @getiposition,...
 'Position',[0.3 0.55 0.15 0.05]);
ivelocityb1=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback', @getivelocity,...
 'Position',[0.3 0.54 0.15 0.05]);
ttimeb1=uicontrol(inpanel,'Style','edit',...
 'String','20','Units','normalized',...
 'BackgroundColor','White','Callback', @gettime,...
 'Position',[0.3 0.5 0.15 0.05]);
align([massb1 stiffnessb1 dampingb1 ipositionb1 ivelocityb1 ...
 ttimeb1],'None','Distribute');

%and sys2
massb2=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback', @getmass2,...
 'Position',[0.83 0.85 0.15 0.05]);
dampingb2=uicontrol(inpanel,'Style','edit',...
 'String','0.2','Units','normalized',...
 'BackgroundColor','White','Callback', @getdamping2,...
 'Position',[0.83 0.57 0.15 0.05]);
stiffnessb2=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback',@getstiffness2,...
 'Position',[0.83 0.56 0.15 0.05]);
ipositionb2=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback', @getiposition2,...
 'Position',[0.83 0.55 0.15 0.05]);
ivelocityb2=uicontrol(inpanel,'Style','edit',...
 'String','1','Units','normalized',...
 'BackgroundColor','White','Callback', @getivelocity2,...
 'Position',[0.83 0.54 0.15 0.05]);
ttimeb2=uicontrol(inpanel,'Style','edit',...
 'String','20','Units','normalized',...
 'BackgroundColor','White','Callback', @gettime2,...
 'Position',[0.83 0.5 0.15 0.05]);
align([massb2 stiffnessb2 dampingb2 ipositionb2 ivelocityb2 ...
 ttimeb2],'None','Distribute');

%Additional parameters for driving force terms sys1
dftypeb=uicontrol(inpanel,'Style','popupmenu',...
 'String',{'None','Sine wave','Pulse'},...
 'Callback', @getdftype,...
 'Units','normalized','BackgroundColor','White',...
 'Value',1,'Position',[0.22 0.4 0.23 0.05]);
dfmagb=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfmag,...
 'Position',[0.3 0.2 0.15 0.05]);
dfperiodb=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfperiod,...
 'Position',[0.3 0.25 0.15 0.05]);
dfsttimeb=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfsttime,...
 'Position',[0.3 0.15 0.15 0.05]);
dfsptimeb=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfsptime,...
 'Position',[0.3 0.1 0.15 0.05]);
align([dftypeb dfmagb dfperiodb dfsttimeb dfsptimeb ],...
 'None','Distribute');

%and sys2
dftypeb2=uicontrol(inpanel,'Style','popupmenu',...
 'String',{'None','Sine wave','Pulse'},...
 'Callback', @getdftype2,...
 'Units','normalized','BackgroundColor','White',...
 'Value',1,'Position',[0.75 0.4 0.23 0.05]);
dfmagb2=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfmag2,...
 'Position',[0.83 0.2 0.15 0.05]);
dfperiodb2=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfperiod2,...
 'Position',[0.83 0.25 0.15 0.05]);
dfsttimeb2=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfsttime2,...
 'Position',[0.83 0.15 0.15 0.05]);
dfsptimeb2=uicontrol(inpanel,'Style','edit',...
 'String','0','Units','normalized',...
 'BackgroundColor','White','Callback', @getdfsptime2,...
 'Position',[0.83 0.1 0.15 0.05]);
align([dftypeb2 dfmagb2 dfperiodb2 dfsttimeb2 dfsptimeb2 ],...
 'None','Distribute');

%Development details
devpanel = uipanel('Parent',fh,'Position',[0.02 0.02 0.22 0.08]);
t(1)={'Developed by and copyright of Dr Martin Gillie'};
t(2)={'University of Manchester'};
t(3)={'Comments to martin.gillie@manchester.ac.uk'};
uicontrol(devpanel,'Style','text',...
 'String',t,'Units','normalized',...
 'HorizontalAlignment','center','Position',[0.02 0.02 0.9 0.9]);
%CALLBACKS FOR GUI. THESE DO THINGS ON EVENTS LIKE TEXT ENTRY IN A
%GUI BOX OR PRESSING A GUI BUTTON

%Function to call solution and plotting when button is pressed
 function solve(hObject,eventdata)
 solution(mass, damping, stiffness, iposition, ivelocity, time, ...
 dftype, dfmag, dfperiod, dfsttime, dfsptime)
 end

% Functions to read values entered in parameter boxes
 function getmass(hObject,eventdata)
 mass(1) = str2double(get(hObject,'string'));
 end
 function getmass2(hObject,eventdata)
 mass(2) = str2double(get(hObject,'string'));
 end
 function getdamping(hObject,eventdata)
 damping{1} = get(hObject,'string');
 damping(1)=strrep(damping(1),'x','y(1)');
 end
 function getdamping2(hObject,eventdata)
 damping{2} = get(hObject,'string');
 damping(2)=strrep(damping(2),'x','y(1)');
 end
 function getstiffness(hObject,eventdata)
 stiffness{1} = get(hObject,'string');
 stiffness(1) = strrep(stiffness(1),'x','y(1)');
 stiffness(1)
 end
 function getstiffness2(hObject,eventdata)
 stiffness{2} = get(hObject,'string');
 stiffness(2) = strrep(stiffness(2),'x','y(1)');
 end
 function gettime(hObject,eventdata)
 time(1) = str2double(get(hObject,'string'));
 end
 function gettime2(hObject,eventdata)
 time(2) = str2double(get(hObject,'string'));
 end
 function getivelocity(hObject,eventdata)
 ivelocity(1) = str2double(get(hObject,'string'));
 end
 function getivelocity2(hObject,eventdata)
 ivelocity(2) = str2double(get(hObject,'string'));
 end
 function getiposition(hObject,eventdata)
 iposition(1) = str2double(get(hObject,'string'));
 end
 function getiposition2(hObject,eventdata)
 iposition(2) = str2double(get(hObject,'string'));
 end
 function getdftype(hObject,eventdata)
 dftype(1) = get(hObject,'Value');
 end
 function getdftype2(hObject,eventdata)
 dftype(2) = get(hObject,'Value');
 end
 function getdfmag(hObject,eventdata)
 dfmag(1) = str2double(get(hObject,'string'));
 end
 function getdfmag2(hObject,eventdata)
 dfmag(2) = str2double(get(hObject,'string'));
 end
 function getdfperiod(hObject,eventdata)
 dfperiod(1) = str2double(get(hObject,'string'));
 end
 function getdfperiod2(hObject,eventdata)
 dfperiod(2) = str2double(get(hObject,'string'));
 end
 function getdfsttime(hObject,eventdata)
 dfsttime(1) = str2double(get(hObject,'string'));
 end
 function getdfsttime2(hObject,eventdata)
 dfsttime(2) = str2double(get(hObject,'string'));
 end
 function getdfsptime(hObject,eventdata)
 dfsptime(1) = str2double(get(hObject,'string'));
 end
 function getdfsptime2(hObject,eventdata)
 dfsptime(2) = str2double(get(hObject,'string'));
 end
end
%%%%%%%%%%%%%%%%%%END OF GUI DEFINITION&amp;&amp;&amp;&amp;&amp;&amp;&amp;&amp;%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%THIS IS THE FUNCTION THAT PLOTS THE DIFFERENTIAL EQUATION SOLUTION

function []=solution(mass,damping,stiffness,iposition,ivelocity,ttime,...
 dftype,dfmag, dfperiod, dfsttime, dfsptime)
%This function solves an ODE for driving force
%initial conditions and time interval. Then it plots various
%graphs in the GUI.

frames=30; %number of frames for plotting
myfunc=@diffeqn;
%Call the ODE solver for two sets of given parameters. The ode itself is
%stored in the function &quot;diffeqn&quot;, handle, myfunc. The calling routine
%is a bit mysterious and copied from an example. Results (y,y',t,f etc)
%are all stored in the output array for later use.

for i=1:2
 [t,y] = ode45(@(t,y) myfunc(t,y,mass(i),damping{i},stiffness{i},...
 dftype(i),dfmag(i), dfperiod(i),dfsttime(i), dfsptime(i)),...
 [0 ttime(i)],[iposition(i) ivelocity(i)]);

 a(i)=size(t,1); %total number of solution increments
 for j=1:a(i)
 output(j,1+4*(i-1))=y(j,1); %position
 output(j,2+4*(i-1))=y(j,2); %velocity
 output(j,3+4*(i-1))=dfforce(t(j), mass(i), dftype(i), dfmag(i),...
 dfperiod(i),dfsttime(i), dfsptime(i));
 output(j,4+4*(i-1))=t(j);
 end
 incs(i)=int16(a(i)/frames); %number of solution increments between frames
end

%set axes dimensions sensibly for plotting later
ysize=zeros(3,2);
for i=1:3
 if max(output(:,i)&gt;=max(output(:,i+4)))
 ysize(i,1)=max(output(:,i))*1.1;
 else
 ysize(i,1)=max(output(:,i+4))*1.1;
 end
 if (min(output(:,i))^2&lt;=(min(output(:,i+4)))^2)
 ysize(i,2)=min(output(:,i+4))*1.1;
 else
 ysize(i,2)=min(output(:,i))*1.1;
 end
 %avoid zero sized axes
 if ysize(i,1)==0
 ysize(i,1)=1;
 end
 if ysize(i,2)==0
 ysize(i,2)=-1;
 end
 %specify balance axes (this can be deleted if desired).
 if ysize(i,1)^2&gt;=ysize(i,2)^2
 ysize(i,2)=-1*ysize(i,1);
 else
 ysize(i,1)=-1*ysize(i,2);
 end
end

if ttime(1)&gt;=ttime(2)
 xsize=ttime(1);
else
 xsize=ttime(2);
end

%Plot information about equations first so that they are immediately
%available. Graphs of solution come below.
for p=1:2
 %setup subplot
 subplot(2,4,4*p,'color',[0 0 0]);
 cla
 plot(0,0); %dummy plot
 axis([-1,1,ysize(1,2),ysize(1,1)],'off');
 if p==1 %select color according to system considered
 col='k';
 else
 col='r';
 end

%Heading for each system considered
 text('units','normalized','position',[0 1.06],'String',...
 ['System ',num2str(p)],'fontsize',10,'fontweight',...
 'bold','color',col);

%GET ALL INFORMATION ABOUT SYSTEM IN FORM OF LATEX STRINGS
 m=num2str(mass(p)); %get mass as latex string

%get driving force terms as latex string
 if dftype(p)==2
 per=num2str(2*pi/dfperiod(p),'%11.3g');
 mag=num2str(dfmag(p),'%11.3g');
 f=[mag,'\sin\left( ',per,'t\right)'];
 elseif dftype(p)==3
 f='\textit{driving force}';
 end
 if (dfsttime(p)~=0 || dfsptime(p)~=ttime(p))
 f='\textit{driving force}';
 end
 if dftype(p)==1
 f=num2str(0);
 end
 if dfsttime(p)==0 &amp;&amp; dfsptime(p)==0
 f=num2str(0);
 end
 %get boundary conditions to pass to latex strings
 ivel=num2str(ivelocity(p));
 ipos=num2str(iposition(p));

%get natural fequency
 if isnan(str2double(stiffness{p})) || isnan(str2double(damping{p}))
 Time=['\textit{period not calculated}'];
 else
 Tp=num2str((2*pi)/(sqrt(4*mass(p)*str2double(stiffness(p))...
 -str2double(damping(p))^2)/(2*mass(p))),'%11.3g');
 Time=['$$T=',Tp,'$$'];
 end

%latex strings for input to text argument. Matlab fussy so can't be
 %put in text command direct
 ddx='\frac{d^2x}{dt^2}';
 dx= '\frac{dx}{dt}';

%Tidy terms before passing to latex fro printing
 damping(p)=strrep(damping(p),'y(1)','x'); %convert back to x as unknown
 stiffness(p)=strrep(stiffness(p),'y(1)','x'); %convert back to x as unknown
 damping(p)=strrep(damping(p),'*',''); %get rid of * symbols
 stiffness(p)=strrep(stiffness(p),'*',''); %get rid of * symbols

 if str2double(damping{p})==0 &amp;&amp; str2double(stiffness{p})~=0
 %miss middle term if no damping
 if str2double(stiffness{p}) &gt; 0 || isnan(str2double(stiffness{p}))
 eqn=['$$',m,ddx,'+',stiffness{p},'x','=',f,'$$'];
 elseif str2double(stiffness{p}) &lt; 0
 eqn=['$$',m,ddx,stiffness{p},'x','=',f,'$$']; %if negative stiffness
 end
 elseif str2double(stiffness{p})==0 &amp;&amp; str2double(damping{p})~=0
 %miss final term if no stiffness
 if str2double(damping{p}) &gt; 0 || isnan(str2double(damping{p}))
 eqn=['$$',m,ddx,'+',damping{p},dx,'=',f,'$$'];
 elseif str2double(damping{p}) &lt; 0
 eqn=['$$',m,ddx,damping{p},dx,'=',f,'$$']; %if negative damping
 end
 elseif str2double(damping{p})==0 &amp;&amp; str2double(stiffness{p})==0
 %or two terms if no stiffness or damping
 eqn=['$$',m,ddx,'=',f,'$$'];
 else
 %include all terms (the normal situation for damped oscillations
 if str2double(stiffness{p}) &gt; 0 &amp;&amp; str2double(damping{p}) &gt; 0
 eqn=['$$',m,ddx,'+',damping{p},dx,'+',stiffness{p},'x','=',f,'$$'];
 elseif isnan(str2double(stiffness{p})) || isnan(str2double(damping{p}))
 eqn=['$$',m,ddx,'+',damping{p},dx,'+',stiffness{p},'x','=',f,'$$'];
 elseif str2double(stiffness{p}) &lt; 0 &amp;&amp; str2double(damping{p}) &lt; 0
 eqn=['$$',m,ddx,damping{p},dx,stiffness{p},'x','=',f,'$$'];
 elseif str2double(stiffness{p}) &gt; 0
 eqn=['$$',m,ddx,damping{p},dx,'+',stiffness{p},'x','=',f,'$$'];
 else
 eqn=['$$',m,ddx,'+',damping{p},dx,stiffness{p},'x','=',f,'$$'];
 end
 end

 %initial conditions as latex string
 icons=['$$','x(0)=',ipos,',\quad ',dx,'(0)=',ivel,'$$'];
 %
 %
 %output latex strings

 text('Interpreter','latex','units','Normalized',...
 'String','System governed by the equation',...
 'Position',[0 0.92],'FontSize',12)
 text('Interpreter','latex','units','Normalized',...
 'String',eqn,'Position',[0 0.73],'FontSize',12)
 text('Interpreter','latex','units','Normalized',...
 'String','with initial conditions',...
 'Position',[0 0.56],'FontSize',12)
 text('Interpreter','latex','units','Normalized',...
 'String',icons,'Position',[0 0.38],'FontSize',12)
 text('Interpreter','latex','units','Normalized',...
 'String','it has a natural period',...
 'Position',[0 0.22],'FontSize',12)
 text('Interpreter','latex','units','Normalized',...
 'String',Time,...
 'Position',[0 0.09],'FontSize',12)
end

%loop through solution for various values of t and take snap shot
%of solution for each t for plotting so graph grows and
%looks like a movie.

timeinc=zeros(2);
for i = 1:frames
 for j=1:2
 timeinc(j)=incs(j)*(i-1)+1; %consider solution up to this increment
 if timeinc(j) &gt; a(j)
 timeinc(j)=a(j); %ensure no overun
 end
 if i==frames
 timeinc(j)=a(j);
 end
 end

 %plots displacement vs time
 subplot(2,4,2);
 plot(output(1:timeinc(2),8),output(1:timeinc(2),5),'color','r');
 hold all
 plot(output(1:timeinc(1),4),output(1:timeinc(1),1),'color','k');
 hold off
 axis([0,xsize,ysize(1,2),ysize(1,1)]);
 line([0 xsize(1)],[0 0],'Color',[0 0 0]);
 ylabel('Position');
 xlabel('Time');
 title('Position against Time','FontWeight','bold');

 %plots velocity vs time
 subplot(2,4,3);
 plot(output(1:timeinc(2),8),output(1:timeinc(2),6),'color','r');
 hold all
 plot(output(1:timeinc(1),4),output(1:timeinc(1),2),'color','k')
 hold off
 axis([0,xsize,ysize(2,2),ysize(2,1)])
 line([0 xsize],[0 0],'Color',[0 0 0]);
 ylabel('Velocity')
 xlabel('Time')
 title('Velocity against Time','FontWeight','bold')

 %Plot driving force
 subplot(2,4,6)
 plot(output(1:timeinc(2),8),output(1:timeinc(2),7),'color','r');
 hold all
 plot(output(1:timeinc(1),4),output(1:timeinc(1),3),'color','k');
 hold off
 line([0 xsize],[0 0],'Color',[0 0 0]);
 xlabel('Time');
 ylabel('Driving force');
 axis([0,xsize,ysize(3,2),ysize(3,1)]);
 title('Driving Force against Time','FontWeight','bold');

 %Plot phase plane
 subplot(2,4,7);
 plot(output(1:timeinc(2),5),output((1:timeinc(2)),6),'color','r');
 hold all
 plot(output(1:timeinc(1),1),output((1:timeinc(1)),2),'color','k');
 hold off
 axis([ysize(1,2),ysize(1,1),ysize(2,2),ysize(2,1)])
 line([ysize(1,2) ysize(1,1) ],[0 0],'Color',[0 0 0]);
 line([0 0],[ysize(2,2) ysize(2,1) ],'Color',[0 0 0]);
 xlabel('Position')
 ylabel('Velocity')
 title('Phase Diagram','FontWeight','bold')
 %pause between frames to ensure behaviour is clear
 pause(0.05)
end

end

%%%%%%%%%%%%%%%%END OF PLOTTING FUNCTION%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%THIS IS THE FUNCTION WHERE THE DIFFERENTIAL EQUATIONS ARE DESCRIBED
%IN MATLAB FORMAT. DIFFERENT EQUATIONS FOR VARIOUS DRIVING FORCE

function [out1] = diffeqn(t,y,mass,damping,stiffness,dftype,...
 dfmag, dfperiod,dfsttime, dfsptime)

f=dfforce(t,mass,dftype, dfmag, dfperiod, dfsttime, dfsptime);
 c=eval(damping); %convert damping from string to number
 k=eval(stiffness);

 out1=zeros(2,1); % a column vector
 out1(1)=y(2); % y(2)=velocity
 out1(2)=(-c/mass)*y(2)-(k/mass)*y(1)+f; % y(1)=position
end

%%%%%%%%%%%%%%%%%FUNCTION EVALUATES DRIVING FORCE AT t%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[f]=dfforce(t,mass,dftype,dfmag, dfperiod, dfsttime, dfsptime)
if dftype==1 %equation with no driving force
 f=0;
elseif dftype==2 %equation with sine wave driving force
 if ((t &gt;= dfsttime) &amp;&amp; (t&lt;=dfsptime))
 f=dfmag*sin(2*pi*t/dfperiod)/mass;
 if dfperiod==0
 f=0; %to avoid division by zero confusing matters
 end
 else
 f=0;
 end
elseif dftype==3 %equation with constant impuse force
 if ((t &gt;= dfsttime) &amp;&amp; (t&lt;=dfsptime))
 f=dfmag/mass;
 else
 f=0;
 end
end
end
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s