Euler's Method and Deflection of Cantilever Beam

Hi, I'm trying to write the code to calculate the deflection of a beam using euler's method. The deflection of the beam is given as w/(24*EI)*(x^4-4Lx^3+6L^2x^2) and the derivative is w/(24*EI)*(4x^3-12Lx^2+12L^2x). This code works except the answers I am looking for are slightly off. So when x = 1.25, y = 3099, but right now I'm getting x = 2.5, y = 3099. Does anyone know why this is or can help me fix this?

w = 21819;
EI = 106;
L = 5; % length
h=1.25; % step's size
N=5; % number of steps
y(1)=1;
for n=1:N       
y(n+1)= y(n)+h*(w/(24*EI))*(4*x(n).^3 - (12*L*x(n).^2) + (12*L^2*x(n)));
x(n+1)=n*h;          % euler's method, using derivative of given function
end
plot(x,y)

Answer 1

[Alan Stevens]

Your derivative is incorrect for your specified deflection equation. I think your deflection equation is also not quite right. Simple Euler is never going to give a good result here because of the rapid variation with x, even with a large number of steps. y(1) must be zero to be consistent with the deflection equation. Have a look at the following:

w = 21819;
EI = 106;
L = 5; % length
N = 100; % number of steps
h = L/N; % step's size
ytrue = @(x)w/(24*EI)*(x.^4 - 4*L*x.^3 + 6*L^2*x.^2); % uniformly loaded cantilever beam
y(1)=0;
x(1) = 0;
for n=1:N-1       
y(n+1)= y(n)+h*w/(24*EI)*(4*x(n).^3 - 12*L*x(n)^.2 + 12*L^2*x(n));
x(n+1)=n*h;          % euler's method, using derivative of given function
end
plot(x,ytrue(x),x,y,'.'),grid
xlabel('x'),ylabel('displacement')
legend('true','Euler')

[JJ]

Hi, sorry I wrote the original equation for the deflection of a beam wrong, so I just fixed that.... The function I wrote in the original code is correct though.. In your code, I'm not getting the answers I confirmed by hand

[Alan Stevens]

You still have y(1) = 1. This isn't consistent with the original deflection equation, where y must equal zero at x = 0. Also note that I used far more steps, so the result at any point will be vastly different from yours.

[JJ]

Yes, I just meant I fixed the error of the wrong equation, not that I edited the code and put it back in the question. I tried the code you suggested and have been playing with it for a while but still reach odd answers. It's very useful, though. thanks.

[Alan Stevens]

You need a hgher order integration routine than simple Euler to get reasonable results here.

[JJ]

It's for an assignment where we just use Euler's method. My point is that the code doesn't match the answers obtained by hand. The problem I am having is that my code results in the correct answers, but for the wrong step. i.e. by hand: when x = 1.25, y = 3099. in Matlab, I'm one step off and the code results in x = 1.25, y = 0, x = 2.5, y = 3099.

When I tried the code you've suggested, it gives completely different answers

[Alan Stevens]

In your hand calculations you almost certainly used x = 1.25 in calculating the right hand side of dydt. This is the backward Euler method and you can incorporate it in your code by changing your loop equations as follows

...
for n=1:N   
x(n+1)=n*h;
y(n+1)= y(n)+h*(w/(24*EI))*(4*x(n+1).^3 - (12*L*x(n+1).^2) + (12*L^2*x(n+1)));
end

[JJ]

Thank you!! All of your responses have been incredibly helpful.

[Ryuxin]

Can we write this program which calculates the N segmenst of beam deflection without using the Euler's method?