%non_e.m
%by Louie L. Yaw - Walla Walla University - 5-8-09
%Example nonlinear explicit analysis of a single dof system, bar in tension
%Demonstration of an incremental EXPLICIT static analysis
%Assume a geometrically nonlinear bar force as f(u)=u+u^3/3
%Therefore k(u)=df/du=1+u^2
%Equal increments are used throughout.
%This is a "load" control type analysis

%Clear all variables
%clear

%%%%%%USER INPUT%%%%%%%%%
%number of increments
ninc=3;
%Max final external load
Pmax=3;
%%%%%%END USER INPUT%%%%%

%%%%%%Explicit Analysis
%The incremental force in each load step
df=Pmax/ninc;
%The current value of the total external applied forces
fg=0;

%The current total global displacement
ug=0;

%Storage arrays for plotting later
ugp=zeros(ninc+1,1);%total global displacements are recorded in this vector
fp=zeros(ninc+1,1);%total global forces are recorded in this vector
%Set first point at the origin(0,0)
ugp(1)=0.0;
fp(1)=0.0;

%Loop over the load increments
for i=1:ninc;
    k=1+ug^2;%calculate the current stiffness, df/du
    du=df/k;%calculate the incremental displacement for load step i
    ug=ug+du;%calculate the current total displacement
    fg=fg+df;%Calculate the current total external force
    ugp(i+1)=ug;%Store the total displacement at the end of increment i
    fp(i+1)=fg;%Store the total external force applied at end of increment i
end;
%%%%%End of Explicit Analysis

%%%%%OUTPUT Final displacements and forces
%The final displacement is
ug
%The total external force is, (or Fext=fg)
Fext=Pmax
%The internal force is
Fint=ug+ug^3/3

%%%%%PLOT exact versus explicit solution%%%%%%
%Plot the explicit and exact force versus displacement solutions
ue=[0:0.1:ceil(max(ugp))]'; %create displacement points for plotting the exact solution
fe=ue+ue.^3./3; %create the array of force data points for plotting exact

%Plot Explicit and exact solutions for comparison
plot(ugp,fp,'bv',ue,fe,'k')
ylabel('Force, F')
xlabel('Displacement, u')
grid on
legend('Explicit','Exact',2)
