%Axial Displacement of Tapered Column.
clear
%Analytical Solution

%Input Variables
Po=1; %Applied load in kips at x=0
P=1;  %Load per foot kips/in along the length of the column
Ao=1;  %Initial area of the column square inches at x=0
L=10; %The length of the column in inches
a=1;  %The parameter 'a' describes how the column cross-section varies linearly
      %based on the formula A(x)=Ao(1+ax/L)
E=29*10^3 ; %The modulus of Elasticity for the column (kips/in^2)

np=50; %The number of data points-1 for x
for i=1:np+1;
    x(i)=(i-1)*L/np;

%The analytical solution of the ODE that describes the column
%The homogeneous solution
uh(i)=(Po*L)/(Ao*E*a)*(log(L+a*x(i))-log(L+a*L));

%The particular solution
xpart(i)=(L+a*x(i))/a*((a*x(i))/(L+a*x(i))*log(L+a*x(i))-log(L+a*x(i)));
lpart(i)=(L+a*L)/a*((a*L)/(L+a*L)*log(L+a*L)-log(L+a*L));
up(i)=(P*L)/(Ao*E*a)*(xpart(i)-lpart(i)+x(i)-L);

%The complete solution is u=uh+up
u(i)=uh(i)+up(i);

%The area at the selected location x(i)
Ax(i)=Ao*(1+a*x(i)/L);

%End of the loop to calculate the area and displacement at every point x
end;

%Calculate stress as a function of x
du=diff(u);
dx=diff(x);
sigmax=du./dx*E;
sigmax(np+1)=sigmax(np)+(sigmax(np)-sigmax(np-1))/dx(np)*abs(dx(np));

%Numerically determined solution goes here
%Specify the number of nodes to be used in the finite difference solution
nnodes=50;
%Define the stiffness matrix k with all zero entries
for i=1:nnodes;
    F(i)=0;
    for j=1:nnodes;
        k(i,j)=0.0;
    end;
end;
%calculate deltax.
deltax=L/(nnodes-1);
%Construct the F and k matrix
for i=1:nnodes;
    %Calculate the value of xi at each node i    
    xi=(i-1)*deltax;
    xn(i)=xi;
    %Calculate some preliminary variables
    lambdai=(L+a*xi)/a;
    %Here we could input P as a function of x
    Pxi=P;
    
    %Construct F and k based on several cases
    %Case 1 if i=1
    if i==1;
        Ri=Pxi*deltax/2;
     %   F(i)=-1.2642;
        F(i)=2*(nnodes-1)*Ri*L/(L*Ao*E)-Po/(Ao*E)*...
            (-2*lambdai/deltax+1);
        k(i,i)=-2*lambdai/deltax^2;
        k(i,i+1)=2*lambdai/deltax^2;
    %Case 3 if i=n-1    
    elseif i == nnodes-1; 
        Ri=Pxi*deltax;
        F(i)=Ri*L*(nnodes-1)/(L*Ao*E);
        k(i,i-1)=lambdai/deltax^2-1/(2*deltax);
        k(i,i)=-2*lambdai/deltax^2;
    %Case 4 if i=n
    elseif i == nnodes;
        k(i,i)=1.0;
        F(i)=0.0;
    %Case 2 if i is greater than 1 and less then n-1    
    else;
        Ri=Pxi*deltax;
        F(i)=Ri*L*(nnodes-1)/(L*Ao*E);
        k(i,i-1)=lambdai/deltax^2-1/(2*deltax);
        k(i,i)=-2*lambdai/deltax^2;
        k(i,i+1)=lambdai/deltax^2+1/(2*deltax);
    end;
end;
%Calculate the numerically determined displacements at each node
un=inv(k)*F';

%Calculate stress as a function of x
ndu=diff(un);
ndx=diff(xn);
nsigmax=ndu./ndx'*E;
nsigmax(nnodes)=nsigmax(nnodes-1)+(nsigmax(nnodes-1)-nsigmax(nnodes-2))...
    /ndx(nnodes-1)*abs(ndx(nnodes-1));
%End of Numerical part

%Plot sigmax vs x
plot(x,sigmax,xn,nsigmax,'k.')

%Graph details
xlabel('Position x (inches)')
ylabel('Axial Stress (ksi)')
title('Axial Stress vs Position For A Tapered Column')
legend('Analytical Solution','Finite Difference Numerical Solution')
grid on

%Pause until the user hits a key
display('please hit a key')
pause

%Plot u vs x
plot(x,-u,xn,-un,'k.')

%Graph details
xlabel('Position x (inches)')
ylabel('Displacement (inches)')
title('Displacement vs Position For A Tapered Column')
legend('Analytical Solution','Finite Difference Numerical Solution')
grid on
%Compare calculated end displacement by virtual work
%to end displacement by the analytical formula above
display('The analytical displacement')
u(1)
display('The numerical displacement')
un(1)
display('The virtual work displacement')
uo=Po*L/(E*Ao*a)*log(1+a)+P*L^2/(E*Ao*a^2)*(a-log(1+a))
display('Numerical Error of Displacement')
ed=abs((u(1)-un(1))/u(1))*100
display('Numerical Error of Stress')
es=abs((sigmax(1)-nsigmax(1))/sigmax(1))*100

