%non_i.m
%by Louie L. Yaw - Walla Walla University - 5-8-09
%Example nonlinear implicit analysis of a single dof system, bar in tension
%Demonstration of an incremental IMPLICIT 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;
%Input Variables for Newton-Raphson iterations
tol=1e-3;%tolerance for equilibrium set by user
maxiter=10;%set maximum number of iterations
%%%%End of User Input

%%%%Start of Implicit 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
ugp(1)=0;
fp(1)=0;

%Loop over the load increments
for i=1:ninc;
    i %print to screen the current load increment
    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
    
    %%%%%Start of statements to include Newton-Raphson iterations
    %Initialize some variables for the Newton-Raphson iterations
    kiter=0;%iteration variable
    dug=0; %Variable of small delta correction to ug to achieve equilibrium
    Fint=ug+ug^3/3; %Calculate the current internal force
    Fext=fg;%Set variable Fext equal to the current external force
    Resid=Fint-Fext%Calculate the residual
    R=sqrt(Resid^2);%Calculate norm of the residual
    
    %Check the residual and max iterations limit, 
    %while above tolerance and max iterations not exceeded do iterations
    while R>tol && kiter<maxiter;
        k=1+(ug+dug)^2; %calculate current stiffness with latest ug+dug
        dug=dug-(1/k)*Resid; %Correct the value of dug based on the residual
        Fint=ug+dug+(ug+dug)^3/3; %Update Fint based on ug+dug
        Resid=Fint-Fext; %Update the residual
        R=sqrt(Resid^2) %Update norm of the residual
        kiter=kiter+1; %Update the number of iterations used to achieve convergence
        %If convergence not achieved and maxiter reached, print a warning
        if kiter==maxiter && R>tol;
            display('WARNING -- Convergence criteria not met')
        end;
    end;
    
    %Output number of iterations required to achieve convergence
    kiter
    ug=ug+dug; 
    %End Newton-Raphson iterations stuff
    
    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 Implicit analysis

%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 the Implicit 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,'bo-',ue,fe,'k')
ylabel('Force, F')
xlabel('Displacement, u')
grid on
legend('Implicit','Exact',2)
