Statically Inderteminant Beam Using Finite Difference Method

Dear all

I am trying to code for the displacements of statically Inderterminant beam using the finite difference method. It appears that I am not that far off from achieving this, however, in one of the points that I have discretized I am getting an incorrect result.

Here, is the graph that comes from the software and what I am trying to replicate:

Here is the graph in MATAB:

As can been seen there is an anomaly at x = 4.1m at which the vertical result is zero. The displacement should only be zero at x = 0m, 5m and 8m.

I don't know if there is an error in my code, but I can't seem to find anything myself.

Here is my code:

clear, clc, close all

% Structural Properties
E = 2.1E+08 % Modulus of elasticity
I = 18890/100^4  % Second moment of area
EI = E*I  % Flexural Rigidity
L = 8  % Beam Length
N = 161 % Number of points

x = linspace(0,L,N)  % Discretise the beam
h = L/(N-1)  % Step size
q = 8 % Load intensity

% Set up matrix using fourth order finite difference

mat1 = diag(ones(1,N)*6) ;                              
mat2 = diag(ones(1,N-1)*-4, 1);
mat3 = diag(ones(1,N-1)*-4, -1);
mat4 = diag(ones(1,N-2)*1, 2);
mat5 = diag(ones(1,N-2)*1, -2);

Mat = [mat1 + mat2 + mat3 + mat4 + mat5]

% Add boundary conditions. Slope = 0 at x = 0 - Use first order finite difference

mat6 = zeros(1,N);
mat6(1,1) = -3;
mat6(1,2) = 4;
mat6(1,3) = -1;
Mat(2,:) = mat6

% Add boundary conditions. Moment = 0 at x = L - using second order finite difference method

mat7 = zeros(1,N);
mat7(1,N-3) = 2
mat7(1,N-2) = -5;
mat7(1,N-1) = 4;
mat7(1,N) = -1
Mat(N-1,:) = mat7

% Remove rows and columns where displacement is zero
% Displacement equals zero at x = 0, x = L and x = 5/8 of beam length

Mat(:,[1,N,(N-1)*5/8]) = []; 
Mat([1,N,(N-1)*5/8],:) = []

% Set up right-hand side

RHS = ones(1,N)*q*h^4/EI;
RHS(2) = 0;
RHS(N-1) = 0

% Remove rows and columns where displacement is zero
% Displacement equals zero at x = 0, x = L and x = 5/8 of beam length

RHS([1,N,(N-1)*5/8]) = []

% Solve for displacement

v = (inv(Mat)*RHS')*1000

% Make new vector with displacements that equal zero

V = [0;v(1:(N+1)/2);0;v(((N+1)/2)+1:end);0]'

% Plot the displacement against beam length

plot(x,-V)
grid on
xlabel('Distance (m)')
ylabel('Displacement (mm)')
title('Local Displacements')

In terms of the code, can anybody see what might be causing this anomaly?

Many thanks in advance!

[Torsten]

V = [0;v(1:(N+1)/2);0;v(((N+1)/2)+1:end);0]'

Obviously, the 0 you insert in the middle is at the wrong position, namely at x(83) = 4.1.

[Scott Banks]

Thanks, Torsten, sorry I didn't see the obvious error.

Only thing now is that the graphs doesn't fall exactly on the 5m mark when the displacement equals zero. Its fall at 4.95m. Is this just an error from the numerical approximation that I will have to live with?

[Torsten]

I guess you "forgot" one discretization point somewhere in between. But since I don't know your mathematical model equations and I didn't dive deep into your discretization, I cannot tell where.

[Scott Banks]

It's okay I have solved it. I needed to use for the displacement of zero at 5m along the beam:

(N-1)*5/8+1