Matlab code: deflection of simply supported beam using ode45

Hi, my code keeps returning a function for a cantilever beam rather than a simply supported. How do I fix this??

My code:

clear; 
clc; 

x0 = 0; 
Y0 = [0; 0; 0; 0]; 
L = 5000;                         % Length of the beam in mm 
xRange = [x0, L]; 

[x, y] = ode45(@(x, y) beamDeflection(x, y, L), xRange, Y0); 
plot(x, y); 
xlabel('Length (mm)'); 
ylabel('Deflection (mm)'); 
title('Simply Supported Beam Deflection'); 

function dydx = beamDeflection(x, y, L) 
E = 69000;                    % Young's modulus of aluminum in MPa 
b = 500;                       % Length of cross-sectional base in mm 
h = 750;                       % Height of cross-sectional height in mm 
I = (b * h^3) / 12;           % Moment of inertia about x 
q = 10000;                     % Load in newtons per mm^2 

if x == 0 || x == L 
    dydx = [0; 0; 0; 0]; 
else 
    dydx = [y(2); y(3); y(4); (q / (E * I))];  
end
end

Answer 1

[Torsten]

Hi, my code keeps returning a function for a cantilever beam rather than a simply supported. How do I fix this??

By using bvp4c instead of ode45. You have a boundary value problem, not an initial value problem because you try to impose boundary condition at x = 0 and x = L.

Here is a symbolic solution:

syms x y(x)
E = 69000;                    % Young's modulus of aluminum in MPa 
b = 500;                       % Length of cross-sectional base in mm 
h = 750;                       % Height of cross-sectional height in mm 
I = (b * h^3) / 12;           % Moment of inertia about x 
q = 10000;                     % Load in newtons per mm^2 
L = 5000;
Dy = diff(y,x);
D2y = diff(y,x,2);
D3y = diff(y,x,3);
D4y = diff(y,x,4);
eqn = D4y == q/(E*L);
conds = [y(0)==0,y(L)==0,D2y(0)==0,D2y(L)==0];
sol = dsolve(eqn,conds)

sol = 

fplot(sol,[0 L])

[Camryn]

Unfortunately it's part of the assessment criteria for me to use ode45. This is still super helpful anyway.

[Torsten]

Then you must estimate y'(0) as y1 and y'''(0) as y3, integrate your system with ode45 with initial conditions [0,y1,0,y3] from x = 0 to x = L, see what comes out at x = L for y(L) and y''(L) (should be both 0), adapt y1 and y3, integrate again, ....

Look up "shooting method" for more details.

[hamawandy]

Why you have not used the moment of inertia in the code ?

[Torsten]

It's a typo. Should be

eqn = D4y == q/(E*I);

instead of

eqn = D4y == q/(E*L);

[hamawandy]

Thank you so much for your help.

I would be so grateful if you let me know how to fill the area under the curve ? which command we can add to do the color fill ?

[Torsten]

syms x y(x)
E = 69000;                    % Young's modulus of aluminum in MPa 
b = 500;                       % Length of cross-sectional base in mm 
h = 750;                       % Height of cross-sectional height in mm 
I = (b * h^3) / 12;           % Moment of inertia about x 
q = 10000;                     % Load in newtons per mm^2 
L = 5000;
Dy = diff(y,x);
D2y = diff(y,x,2);
D3y = diff(y,x,3);
D4y = diff(y,x,4);
eqn = D4y == q/(E*I);
conds = [y(0)==0,y(L)==0,D2y(0)==0,D2y(L)==0];
sol = dsolve(eqn,conds);
xn = 0:0.1:L;
yn = subs(sol,x,xn);
area(xn,yn)

●

[hamawandy]

Thank you so much ...... I am so grateful to you.

Answer 2

[hamawandy]

Dear Sir

I am facing a problem in sketching the shear force diagram from the span 4m till 8m.......the line should be horizontal from 4m to 8m.

Can you help me to fix this error ?

This is the code:

R1=6;
R2=2;
X= linspace (0, 8, 100);
V= (R1-2*X).*(X>=0) - (R2)*(X>=4);
subplot (211)
plot (X, V)
area (X,V, 'FaceColor', 'r', 'LineWidth', 2);

[Torsten]

Maybe

V= (R1-2*X).*(X>=0).*(X<=4) - (R2)*(X>=4);

or

V= (R1-2*X).*(X>=0).*(X<4) - (R2)*(X>=4);

instead of

V= (R1-2*X).*(X>=0) - (R2)*(X>=4);

?

[hamawandy]

Wow......it is done

Finally, what about the bending moment diagram equation ? should I multiply by the distance for each term ? or need a new code ?

[Torsten]

Which function for M do you want to plot in a mathematical notation (not code) ?

[hamawandy]

Thank you so much I have fixed the error just now.

The final code for both shear force and bending moment diagrams is the following:

R1=6;
R2=2;
X= linspace (0, 8, 100);
V= (R1-2*X).*(X>=0).*(X<=4) - (R2)*(X>=4);
subplot (211) 
plot (X, V)
area (X,V, 'FaceColor', 'r', 'LineWidth', 2);
M=(R1*X-X.^2 ).*(X>=0).*(X<4) + (R1*X-8*(X-2)).* (X>=4).* (X<=8);
subplot(212)
plot (X,M)
area(X,M, 'FaceColor', 'b', 'LineWidth', 2);
xline(0, 'Color', 'k', 'LineWidth', 2); % Draw line for Y axis
yline(0, 'Color', 'k', 'LineWidth', 2); % Draw line for X axis

●