12 views (last 30 days)

Hey guys, I am unable to make my program get a result of w as a function of EI and L

syms w x EI A B C D L

x0=0;

eqns = [EI*diff(w,3)==heaviside(x-x0)+A, EI*diff(w,2)==heaviside(x-x0)*(x-x0)+A*x+B, EI*diff(w)==heaviside(x-x0)*(x-x0)^2/2+A*x^2/2+B*x+C, EI*w==heaviside(x-x0)*(x-x0)^3/6+A*x^3/6+B*x^2/2+C*x+D, w(x==L/2)==0, diff(w(x==-L/2))==0];

ctes = [A B C D];

[A, B, C, D]=solve(eqns, ctes);

w(x==0)

This is what I've coded. As you can see it is 4 equations and 2 boundary conditions (The 2 last ones) I want it to solve the 4 constants A, B, C and D, so I get the full result of the function w(x) with those Constants solved and then finally evaluate w for x=0 (in this case). But it is giving me this result:

ans =

[ empty sym ]

I am pretty sure it is something about syntax in the equations definition, but I can't find it. Could you give it a look guys? I would really appreciate. This example should give the following solution w(x=0) = L^3/(192*EI)

Thank you very much!

John D'Errico
on 27 Jun 2020

These are differential equations, a boundary value problem.

We use dsolve to solve it. For example, a simple example is:

syms x(t)

Eq = diff(x,t,2) == 1 + t;

dsolve(Eq,x(0) == 0,x(1) == 2)

ans =

(t*(t^2 + 3*t + 8))/6

Now, what happens when you just throw it at solve?

solve(Eq,x(0) == 0,x(1) == 2)

Warning: Unable to find explicit solution. For options, see help.

> In sym/solve (line 317)

ans =

Empty sym: 0-by-1

Finally, even though you used solve to try to find a solution, that last line, w(x==0) does not do what you think it does. In fact, in your code, even though you tried to return A,B,C,D, those values are not magically inserted into the relationship for w.

Oh, this is not how you create a boundary condition:

w(x==L/2)==0

Again, that does NOT set x to a value of L/2, and then substitute it into w. Instead, you would formulate that boundary condition for dsolve as

w(L/2) == 0

So, first, solve the ODE as a differential equation. Use dsolve.

John D'Errico
on 27 Jun 2020

While I know the wiki page claims the two ways of formulating the problem are equivalent,

it seems possible that dsolve is not happy with using a dirac delta at an end point to formulate the problem. For example, reading the answer to this question:

What you tried

syms w(x) x EI L keq

x0=L;

assume(L,'positive')

EDO = EI*diff(w,x,4) == dirac(x-x0);

Dw = diff(w,x);

DDw = diff(w,x,2);

DDDw = diff(w,x,3);

CC = [w(0) == 0, Dw(0) == 0, DDw(L) == 0, DDDw(L) == 0];

must fail though. As even the euler-bernouli page suggests it will. See this statement on the wiki page:

"Note that shear force boundary condition (third derivative) is removed, otherwise there would be a contradiction."

I'll admit that I tried a few things with no success to get dsolve to happily use dirac to implement this as an impulse response at the right hand end point.

So here, with 3 boundary conditions. Strangly, I could not get dsolve to solve this over the interval [0,L]. So here I solved it over a unit interval, and then transformed the problem. (I was also lazy and implicitly left EI == 1. That should not be very significant, except to make things a bit simpler to read.)

W1 = dsolve(diff(w,4) == dirac(x - 1),w(0) == 0,dw(0) == 0,ddw(1) == 0)

W1 =

x*(sign(x - 1)/4 + 1/4) - sign(x - 1)/12 - x^2*(sign(x - 1)/4 - C1 + 1/4) + x^3*(sign(x - 1)/12 - C1/3 + 1/12) - 1/12

As you should see, there is a constant of integration remaining in there: C1.

Now, transform to the interval [0,L].

syms u

Wu = subs(W1,x,u/L)

Wu =

(u*(sign(u/L - 1)/4 + 1/4))/L - sign(u/L - 1)/12 - (u^2*(sign(u/L - 1)/4 - C1 + 1/4))/L^2 + (u^3*(sign(u/L - 1)/12 - C1/3 + 1/12))/L^3 - 1/12

The sign function in there can be viewed as a heaviside function.

Does it satisfy the boundary conditions?

>> subs(Wu,u,0)

ans =

0

>> subs(diff(Wu,1),u,0)

ans =

0

>> subs(diff(Wu,2),u,L)

ans =

NaN

>> limit(diff(Wu,2),u,L,'left')

ans =

0

As you can see, on the right end I needed to take the limit as u --> L from the left.

If I now take the limit of the third derivative, as u approaches L from the left side...

>> limit(diff(Wu,3),u,L,'left')

ans =

-(2*C1)/L^3

I get something that is consistent with the use of an impulse response on that left end.

But honestly, I felt I had to work much harder to get to this point than I did by just imposing a shear load directly in the boundary conditions. I'm not sure it is worth it. :)

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.