Create a Script to solve for beam deflection

% March 17

%% Example 1, use ode45 to solve for a system of 3 ODEs

clear,clc

close all

[t,y] = ode45(@dydtfun,[0,12],[0,1,1]);

plot(t,y(:,1),t,y(:,2),t,y(:,3))

xlabel('time')

ylabel('y values')

legend('y1','y2','y3')

%% Example 2, use ode45 to solve for a vibration system

y0 = [0.5,0];

tspan = [0,50];

[t,y] = ode45(@vibfun,tspan,y0);

figure

subplot(2,1,1)

plot(t,y(:,1))

subplot(2,1,2)

plot(t,y(:,2))

%% Example 2 modified, compare ode45 and ode15s

clear,clc

close all

y0 = [0.5, 0];

tspan = [0,1e6]; % 1e6 is 10^6

tic

[t,y] = ode45(@vibsfun,tspan,y0);

toc

tic

[t1,y1] = ode15s(@vibsfun,tspan,y0);

toc

plot(t,y(:,1),t,y(:,2),t1,y1(:,1),'o',t1,y1(:,2),'o')

legend('x by ode45','v by ode45','x by ode15s','v by ode15s')

%% Define all the ode functions

% Example 1

function yprime = dydtfun(t,y)

yprime = [0;0;0];

yprime(1) = y(2)*y(3)*t;

yprime(2) = -y(1)*y(3);

yprime(3) = -0.51*y(1)*y(2);

end

% alternatively

% yprime = [y(2)*y(3)*t;-y(1)*y(3);-0.51*y(1)*y(2)]

% Example 2

function yprime = vibfun(t,y)

% y(1) is x

% y(2) is xdot

% yprime(1) is xdot or y(2)

% yprime(2) is x2dot

m = 10;

k = 4;

c = 2;

f0 = 0.05;

w = 2;

yprime = [0;0];

yprime(1) = y(2);

yprime(2) = f0/m*sin(w*t) - c/m*y(2) - k/m*y(1);

end

% Example 2 modified to be a stiff ode problem

function yprime = vibsfun(t,y)

% y(1) is x

% y(2) is xdot

% yprime(1) is xdot or y(2)

% yprime(2) is x2dot

m = 1;

k = 1e-3;

c = 1;

yprime = [0;0];

yprime(1) = y(2);

yprime(2) = - c/m*y(2) - k/m*y(1);

end

% Mar. 22

%% Beam deflection

% define all the known parameters

% divide the length into N subdivisions

%

N = 10000;

dx = L/N;

N1 = round(L1/dx);

N2 = round((L2-L1)/dx);

% calculate the reactions forces

R1 =

R2 =

% create the known vector, p_n

Pn = zeros(N+1,1);

% replace the elements in Pn using the specific expression/values for that

% beam section

for n = 2:N1+1

Pn(n) =

end

for n = N1+2:N1+N2+1

Pn(n) =

end

for n = N1+N2+2:N

Pn(n) =

end

% define the coefficient matrix A

A = zeros(N+1,N+1);

% use diag function to modify the matrix A, matching the tridiagnol matrix

% solve for the linear equations

y = A\Pn;

Answer 1

[Himanshu]

Hey Deric,

I understand that you are trying to write a MATLAB script to solve for beam deflection. Please go through the following documentation page to find examples related to the same:

https://www.mathworks.com/matlabcentral/fileexchange/113670-beam-bending-and-deflection

The following MATLAB Answers thread might also be useful:

• Euler's Method and Deflection of Cantilever Beam- https://www.mathworks.com/matlabcentral/answers/648468-euler-s-method-and-deflection-of-cantilever-beam

Hope this helps!