Solving 4th order ode

7 visualizaciones (últimos 30 días)
marijn de volder
marijn de volder el 15 de Jul. de 2016
Comentada: Star Strider el 16 de Jul. de 2016
I try to solve the ODE shown on the figure. E and I are constants and q is a point Load in the middle of the beam i.e. at x=5. I want to obtain the settlement (=s) in function of x. The beam starts at x=0 and ends at x=10.
  • My first problem is to write the correct fuction for q. I used: q*heaviside(x-5)
  • The second problem is solving the ODE. I used : dsolve function. I have 4 Boundary conditions:
  1. The Moment is 0 at x=0
  2. The Moment is 0 at x=10
  3. The Shear force is 0 at x=0
  4. The Shear force is 0 at x=10
This is my code. Thanks in advance and keep in mind that i am a beginner in matlab.
EI=10;
k=4;
q=1;
x=0:0.1:10;
syms s(x);
Ds = diff(s);
D2s = diff(s,2);
D3s = diff(s,3);
dsolve(q*heaviside(x-5)==diff(s,4)*EI+k*s,D2s(0)*EI==0,D2s(10)*EI==0,D3s(0)*EI==0,D3s(10)*EI==0)

Respuesta aceptada

Star Strider
Star Strider el 15 de Jul. de 2016
Mechanical engineering is nto my area of expertise, so I won’t comment with respect to your using heaviside correctly here.
You seem to have set it up correctly, but it’s not giving the output in the image you posted:
syms s(x);
EI=sym(10);
k=sym(4);
q=sym(1);
Ds = diff(s);
D2s = diff(s,2);
D3s = diff(s,3);
s = dsolve(q*heaviside(x-5)==diff(s,4)*EI+k*s,D2s(0)*EI==0,D2s(10)*EI==0,D3s(0)*EI==0,D3s(10)*EI==0);
s = simplify(s, 'steps',10)
figure(1)
subplot(4,1,1)
ezplot(imag(s), [0 10])
title('s(x)')
subplot(4,1,2)
ezplot(imag(diff(s)), [0 10])
title('\theta')
subplot(4,1,3)
ezplot(EI*imag(diff(s,2)), [0 10])
title('M')
subplot(4,1,4)
ezplot(EI*imag(diff(s,3)), [0 10])
title('T')
It simplifies to:
s =
- (sign(x - 5)*1i)/8 - 1i/8
You can take out the simplify call if you want, but it doesn’t significantly affect the results.
  2 comentarios
marijn de volder
marijn de volder el 16 de Jul. de 2016
Thanks a lot, this really helped me. I also saw that my function for the force should be: q*( heaviside(x-4.999)-heaviside(x-5.001)). Unfortunately I still don't get the correct result, but i'am getting closer.
Star Strider
Star Strider el 16 de Jul. de 2016
My pleasure.
I’m certain you will get it sorted. Unfortunately, I can’t provide any insight other than with respect to the code.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by