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!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! 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! matrices, finite diffference method, boundary value problems MATLAB Answers — New Questions

​