Several errors in my code.

Hi there,

We were working on this code in class and I have been informed that there are several errors in my code below, but for the life of me I can't figure out where!

function res = FEM_1D(L, k, f, xf, S, N)
% L is vector of section length
% k is vector of stiffness AE/L
% f is vector of force magnitude
% xf is vector of force location
% S is vector of supports
% N is vector of elements in a section (optional)
if nargin <6 %check if the optional parameter N is not given
    N = 10*ones(1,length(L)); % Make a matrix (10,10,10...10)
end
n = dot(L,N)+1; % Number of nodes
X = zeros(n,1); %Number of node coordinates
F = zeros(n,1); %Vector force on each node
res.U = zeros(n,1); %Vector of displacements, but stored in res
K = zeros(n,n); %Global stiffness matrix
count=1; %my counter
for i =1:length(L) %Loop over all sections of the 1D bar
    for j=1:N(i) %Loop over all elements in each section
        count=count + 1; %counting
        X(count,1)= X(count-1,1)+L(i)/N(i); %coordinates
        ke=k(i)*N(i)*[1 -1;-1 1]; %Local matrix of stiffness
        K(count-1:count,count-1:count)=K(count-1:count,count-1:count)+ke;
        %add local matrix to the global matrix by superposition
        %K(count-1:count,count-1:count)=K(count-1:count,count-1:count)+ke(1,1);
        %K(count-1:count,count-1:count)=K(count-1:count,count-1:count)+ke(1,2);
        %K(count-1:count,count-1:count)=K(count-1:count,count-1:count)+ke(2,1);
        %K(count-1:count,count-1:count)=K(count-1:count,count-1:count)+ke(2,2)
    end
end
forceid = interp1(X, 1:length(X),xf,'nearest','extrap');
%find the coordinate of nodes that is closest to the force location
F(forceid)=f; %apply force for the nodes found
supportid = interp1(X,1:length(X),S,'nearest','extrap'); 
%find the coordinate of nodes that is closest to the support location
movable = setdiff(1:n,supportid);
Kreduced = K(movable,movable); %copy only movable value
Freduced = F(movable); %copy only the movable force value
Ureduced = Kreduced\Freduced; %calculate the displacement U= F/k
res.U(movable)=Ureduced; %export value to res
end

It runs the script fine and I believe the answer it spits out are correct, but would there be any suggestions on how to improve?

Any help would be greatly appreciated!

[Bhaskar R]

Can you provide inputs (L, k, f, xf, S, N)?

[Thomas Stokes]

Sorry about that:

L=1, k=10000, f=100, xf=0.5, s=0.2 and N is an unknown.

[Mario Malic]

Check if your stiffness matrix is symetric.

I don't understand this "reduced" something. Once you get the stiffness matrix, invert it and multiply with multiply with force vector to get the displacements.