Desperately seeking help with solving a differential equation

[hushish]

Hi,

I desperately need help to solve the following differential equation for buckling of a beam with a uniform axially applied force and a point force:

∂y(x)²/∂x²+(P+Q.x).y(x)=0
Where P and Q are constants. P is known and Q is the critical axial uniform force (N/mm) that will cause buckling.

Regards and thanks,

Hushish

 

[JJacquelin]

The change of variable X=(P+Q.x).Q^(-2/3) leads to y''(X)=X.y which is an Airy ODE.
Solutions are : y=c1.Ai(X)+c2.Bi(X)
Ai and Bi are the Airy functions.

 

[hushish]

Thanks JJacquelin.

How did you change the variable? Please post steps and excuse my maths ignorance.

Would Ai(X) and Bi(X) be simply the standard Airy functions with (P+Q.x).Q^(-2/3) substituted for X?

 

[Bill Simpson]

This will not show you the steps, but you can compare his fine result with this
http://www.wolframalpha.com/input/?i=DSolve%5By%27%27%5Bx%5D+%2B+%28p+%2B+q*x%29+y%5Bx%5D+%3D%3D+0%2C+y%5Bx%5D%2C+x%5D
Then you can compare this
http://www.wolframalpha.com/input/?i=DSolve%5By%27%27%5Bx%5D+%2B+x+y%5Bx%5D+%3D%3D+0%2C+y%5Bx%5D%2C+x%5D
with this
http://www.wolframalpha.com/input/?i=DSolve%5By%27%27%5Bx%5D+%2B+q+x+y%5Bx%5D+%3D%3D+0%2C+y%5Bx%5D%2C+x%5D
to try to guess how he was introducing the q into the denominator to get your solution.

Then you can think about how p plays into this.

You can check their definition of Airyai and Airybi here
http://reference.wolfram.com/mathematica/ref/AiryAi.html
http://reference.wolfram.com/mathematica/ref/AiryBi.html
to make certain that they are not using some different definition of those. (Click on "Details" on each of those to see the information needed)

 

[hushish]

Bill,

Thanks for the info. I plugged the equation into Matlab and got the following solution for y''+q.x.y=0:

y =C1*AiryAi(-(q/E/I)^(1/3)*x)+C2*AiryBi(-(q/E/I)^(1/3)*x)

From this, I input (-(q/E/I)^(1/3)*x) into the Airy's functions. y(0)=0, therefore C1=0, and y(L)=0, therefore I solve for AiryBi = 0 and obtain q.

I have 2 further questions:
1) If I want to obtain y''(L)=0, and I substitute the values into the Airy function for y(x) and differentiate twice, I don't obtain what I should get. In Matlab, y''(x)=C1*AiryAi(-(q/E/I)^(1/3)*x)+C2*AiryBi(-(q/E/I)^(1/3)*x), which is the same function as for y(x). Is that correct?

2) y''+(P+q.x).y=0 yields the folliowing solution in Matlab:
y(x)=C1*AiryAi(-(q/E/I)^(1/3)*(P+q*x)/q)+C2*AiryBi(-(q/E/I)^(1/3)*(P+q*x)/q)
However, by varying P to obtain q in the same way as illustarted above, I get close, but I never obtain the solution y(x)=0.
Am I doing something wrong?

Thanks again,

 

[Bill Simpson]

Bill,

Thanks for the info. I plugged the equation into Matlab and got the following solution for y''+q.x.y=0:

y =C1*AiryAi(-(q/E/I)^(1/3)*x)+C2*AiryBi(-(q/E/I)^(1/3)*x)

I get what seems to be about the same result from Mathematica, but I am not sure where your E and I are from.

From this, I input (-(q/E/I)^(1/3)*x) into the Airy's functions. y(0)=0, therefore C1=0, and y(L)=0, therefore I solve for AiryBi = 0 and obtain q.

I usually give initial conditions with the differential equation and let the tool solve for the constants. I seem to make fewer mistakes when I do that.

When I include the condition y(0)=0 that eliminates one of the constants and I get

y[x]= AiryAi[(-q)^(1/3) x] C[1] - (AiryBi[(-q)^(1/3) x] C[1])/Sqrt[3]

but that does not seem to give the same answer that you have found.

When I try to include both y(0)=0 and y(L)=0 I get the solution y(x)=0.

I have 2 further questions:
1) If I want to obtain y''(L)=0, and I substitute the values into the Airy function for y(x) and differentiate twice, I don't obtain what I should get. In Matlab, y''(x)=C1*AiryAi(-(q/E/I)^(1/3)*x)+C2*AiryBi(-(q/E/I)^(1/3)*x), which is the same function as for y(x). Is that correct?

I think I am getting lost here because I can't reproduce your y(L) and so can't guess what y''(L) is.

2) y''+(P+q.x).y=0 yields the folliowing solution in Matlab:
y(x)=C1*AiryAi(-(q/E/I)^(1/3)*(P+q*x)/q)+C2*AiryBi(-(q/E/I)^(1/3)*(P+q*x)/q)
However, by varying P to obtain q in the same way as illustarted above, I get close, but I never obtain the solution y(x)=0.
Am I doing something wrong?

I have no idea. I think we have to show every step in detail and try to figure out why you are getting something very different from what I am getting. Maybe then we can try to figure out whether your result just looks different or whether one or the other is an error.

Thanks again,

Glad I could help

In[1]:= Simplify[DSolve[{y''[x] + q x y[x] == 0}, y[x], x]]

Out[1]= {{y[x] -> AiryAi[(-q)^(1/3) x] C[1] + AiryBi[(-q)^(1/3) x] C[2]}}

In[2]:= Simplify[DSolve[{y''[x] + q x y[x] == 0, y[0] == 0}, y[x], x]]

Out[2]= {{y[x] -> AiryAi[(-q)^(1/3) x] C[1] - (AiryBi[(-q)^(1/3) x] C[1])/Sqrt[3]}}

Check the result is 0 when x == 0

In[3]:= ReplaceAll[AiryAi[(-q)^(1/3) x] C[1]-(AiryBi[(-q)^(1/3) x] C[1])/Sqrt[3],x->0]

Out[3]= 0

In[4]:= Simplify[DSolve[{y''[x] + q x y[x] == 0, y[0] == 0, y[L] == 0}, y[x], x]]

Out[4]= {{y[x] -> 0}}

 

[hushish]

Hi Bill,

Sorry for the inclusion of the E and I term without explanation. The problem is for buckling of a simply supported beam with an intital force P and a distributed load q. Now the answer converges when there is no P term. However, if there is an initial P term, I am left with the following solution:

y =C1*AiryAi((-P-q*x)/q*(q/E/I)^(1/3))+C2*AiryBi((-P-q*x)/q*(q/E/I)^(1/3))

Since y(0)=0, C1=0. y(L) = 0, therefore C2*AiryBi((-P-q*L)/q*(q/E/I)^(1/3))=0

If I expand the terms:

y(x)=C2*(X+X^4/12+X^7/504+X^10/45360+...)

Where:

X=((-P-q*x)/q*(q/E/I)^(1/3))

The problem is that for any real value of P and q, y(0)=0 does not converge. This is the problem that I am facing.

Thanks again,

hushish

 

[Bill Simpson]

However, if there is an initial P term, I am left with the following solution:

y =C1*AiryAi((-P-q*x)/q*(q/E/I)^(1/3))+C2*AiryBi((-P-q*x)/q*(q/E/I)^(1/3))

Since y(0)=0, C1=0

I don't think I find that C1=0 when y(0)=0.

In[1]:= Simplify[DSolve[{y''[x]+(p+q x)y[x] == 0}, y[x], x]]

Out[1]= {{y[x]->AiryAi[-((p+q x)/(-q)^(2/3))]C[1]+AiryBi[-((p+q x)/(-q)^(2/3))]C[2]}}

Now introduce y(0)=0

In[2]:= Simplify[DSolve[{y''[x]+(p+q x)y[x] == 0, y[0] == 0}, y[x], x]]

Out[2]= {{y[x] -> (AiryAi[-((p + q x)/(-q)^(2/3))] -
      (AiryAi[(p q)/(-q)^(5/3)] AiryBi[-((p + q x)/(-q)^(2/3))])/
      AiryBi[(p q)/(-q)^(5/3)]) C[1]}}

Which is sort of like saying c1=1 and
c2=AiryAi[(p q)/(-q)^(5/3)]/AiryBi[(p q)/(-q)^(5/3)] and
then any scalar multiple of that solution is still a solution.

Now check that really is 0 when x=0

In[3]:= Simplify[ReplaceAll[Out[2], x -> 0]]

Out[3]= 0

Maybe you can find a way to show that is just AiryBi[-((p + q x)/(-q)^(2/3))] C[2] but I can't seem to do that. So I suspect there may be a small flaw in your derivation that c1=0.

Can you check that?

Can you try giving initial conditions to Matlab and let it calculate the constants?

If I then add the second condition that y(L)=0 I get only

In[4]:= Simplify[DSolve[{y''[x]+(p+q x)y[x] == 0, y[0] == 0, y[L] == 0}, y[x], x]]

Out[4]= {{y[x] -> 0}}

which might be a hint why it doesn't converge and you cannot find a solution, because there isn't one as the problem is formulated. Note: This doesn't change if I replace L with 1 to try to make it easier.

EDIT

If I pick values like p=5 and q=3 and plot the real and complex components of the solution I found above for y(0)=0 then there are a sequence of points where re=im=0. In other words, it isn't that there are no solutions for any x other than 0, there are some. Thus if your L happened to be exactly one of those values then there should be a solution. The first three roots for p=5 and q=3 are at about 1.20258, 2.19111, 3.06609. All this hinges on the handling of roots of negative numbers. And I don't know if it really applies to how your beam really behaves.

EDIT

Telling Mathematica that p,q,x,L are all >=0 allows it to simplify slightly more

In[2]:= Assuming[p > 0 && q > 0 && x >= 0, 
 Simplify[DSolve[{y''[x] + (p + q x) y[x] == 0, y[0] == 0}, y[x], x]]
 ]

Out[2]= {{y[x] -> (AiryAi[((-1)^(1/3) (p + q x))/q^(2/3)] -
      (AiryAi[(p q)/(-q)^(5/3)] AiryBi[((-1)^(1/3)(p+q x))/q^(2/3)])/AiryBi[(p q)/(-q)^(5/3)]) C[1]}}

Using the default cube root of -1 gives the same result and plot as before.
Forcing the cube root to be -1 gives a complex expression never equal to zero for any x<>0.

 

[hushish]

Hi Bill,

You are right, C1 does not equal zero. Nonetheless my answers still do not converge to the correct results. Since the buckling equation is a transcendental equation, it cannot be solved explicitly. Nevertheless, the values of q that satisfy the equation can be determined numerically. I try to isolate the problem with boundary conditions so that C1 is given in terms of C2, and then solve inside the brackets by iteration to obtain q.

Just out of interest, if you input the following values into Mathematica, can you get the following value for q:

y''[x]+1/(E*I)*(p+q x)y[x] = 0
y[0]=0
y[L]=0
E=1000
I=1000
L=100
P=500
q should be around 10.

Thanks,

hushish

 

[Bill Simpson]

In[1]:= e = 1000;
i = 1000;
L = 100;
DSolve[{500 y''[x] + 1/(e*i)*(p + q x) y[x] == 0, y[0] == 0, y[L] == 0}, y[x], x]

Out[4]= {{y[x] -> 0}}

As I've said a few times, it is firmly convinced there is no general solution other than y[x]=0, but some of the decisions they made decades ago means that it can miss isolated points that satisfy when it is being asked for a general solution.

Can you ask the question again in a different way and maybe I can find a way to get it to look for your 10 when p is unknown.

 

[hushish]

Hi Bill,

The trivial solution will always tend to C1=C2=0 for buckling, as C1 and C2 are irrelevant. It's always the terms inside the bracket that determine bthe critical buckling load. Only once I have an expression for y[x] can I manipulate the boundary conditions to obtain q.

I am trying to recreate the attached graph from Michael Niu, Airframe Stress Analysis and Sizing.

Thanks again,

 

[Bill Simpson]

I think you have reached the limits of my skills on this. Are there any other places you could ask about this?

 

[hushish]

Thanks Bill for your time and help;, I will post something here if I find anything.

 

[hushish]

Hi Bill,

For some reason, the approach I took (with your help) above seems to work and my answers are now converging.

Perhaps you can add some insight into the next problem:

If I adjust the boundary conditions, the problem now requires a particular solution:

y''[x]/(E*I)+(P+q x)*y[x] = R*(L-x)

Where R is the reaction at the boundary conditions, and L is the beam's length. If I put this query into Wolfram I get a huge solution. In Matlab I get the following:

y(x)=
AiryAi((-P-q*x)/q*(q*E*I)^(1/3))*C2+AiryBi((-P-q*x)/q*(q*E*I)^(1/3))*C1+R*E*I*pi*(-Int(AiryBi(-(q*E*I)^(1/3)*(P+q*x)/q)*(-L+x),x)*AiryAi(-(q*E*I)^(1/3)*(P+q*x)/q)+Int(AiryAi(-(q*E*I)^(1/3)*(P+q*x)/q)*(-L+x),x)*AiryBi(-(q*E*I)^(1/3)*(P+q*x)/q))/(q*E*I)^(1/3)

My question is this:

Can I inegrate an Airy function Int(AiryAi(-(q*E*I)^(1/3)*(P+q*x)/q)*(-L+x),x) once I have it in the series solution form, or is it more complex than that?

Regards,

 

[Bill Simpson]

I think you are falling into a very common trap, you may be assuming that I have been hiding in the room behind you and I am aware not only of every keystroke, but of all the thoughts inside your head.

Since I am actually not aware of any of that, can you perhaps give me a hint of which of your variables are complex? Which are real? Which are greater than zero? What are these "boundary conditions" that you have "adjusted" and how have you adjusted them? Are any of these things numbers that might be revealed? Or does all this have to be done strictly symbolically with no domain information that can be leaked to the outside? etc.

This

Simplify[DSolve[y''[x]/(e*i) + (P + q x)*y[x] == R*(L - x), y[x], x]]

seems to give a fairly simple expression with no conditions and all in terms of your unknown coefficients. To do better than that I, and certainly Mathematica, need more specifics.

Just making a wild guess and telling it

Assuming[e > 0 && i > 0 && P > 0 && q > 0 && r > 0 && L > 0 && x > 0, 
 Simplify[DSolve[y''[x]/(e*i) + (P + q x)*y[x] == R*(L - x), y[x], x]]]

doesn't seem to help at all.

 

[hushish]

Hi Bill,

Sorry about that. The boundary conditions have changed in that now rotations and deflections are zero at x=0. Thus y(0)=0 and y'(0)=0. y(L)=0. R is the reaction force at x=0 and is real and not a function of x. All the assumptions you have expressed in Mathematica are correct. I am just looking for the particular solution to the R*(L-x) term so that I can have an expression for y(x). In Matlab, it's not as simple as R/P*(L-x). Then by applying boundary conditions, I can eliminate R and obtain the critical buckling load q*L. This procedure is highlighted in the attached example.



Thanks,

hushish

 

[Bill Simpson]

Ah, finally the information needed to be able to figure all of this out.

In[1]:= DSolve[v''[x] + k^2 v[x] == R/(e i) (L - x), v[x], x]

Out[1]= {{v[x] -> (L R - R x)/(e i k^2) + C[1] Cos[k x] + C[2] Sin[k x]}}

Manually substitute k^2 e i == P because it is easier to do that manually than it is to fight with Mathematica

C[1] Cos[k x] + C[2] Sin[k x] + R/P (L - x)

That matches the equation 11-39 in your text.

Note: Mathematica assigns C[1] and C[2] in the opposite order from the text.

Now see whether automating the finding of the first coefficients works, because trying to apply all three boundary conditions always fails

In[2]:= Apart[DSolve[{v''[x] + k^2 v[x]==R/(e i) (L-x), v[0]==0, v'[0]==0},v[x],x]]

Out[2]= {{v[x] -> -((R (-L + x + L Cos[k x]))/(e i k^2)) + (R Sin[k x])/(e i k^3)}}

Which is

R/P (L - x) - (R L)/P Cos[k x] + R/P Sin[k x]/k

I think I see the C1 and C2 in that, but trying to solve for R when x==L gives

R/P (L - L) - (R L)/P Cos[k x] + R/P Sin[k x]/k==0

And I suppose that gives

(R L)/P Cos[k x] + R/P Sin[k x]/k==0

R(L Cos[k x] + Sin[k x]/k)==0

But I don't understand what he expects to do with that to find R.

Now try doing those manually.

v[0]==0

In[3]:= (C[1] Cos[k x] + C[2] Sin[k x] + R/P (L - x) /. x -> 0) == 0

Out[3]= (L R)/P + C[1] == 0

v'[0]==0

In[4]:= (D[C[1] Cos[k x] + C[2] Sin[k x] + R/P (L - x), x] /. x -> 0) == 0

Out[4]= -(R/P) + k C[2] == 0

v[L]==0

In[5]:= (C[1] Cos[k x] + C[2] Sin[k x] + R/P (L - x) /. x -> L) == 0

Out[5]= C[1] Cos[k L] + C[2] Sin[k L] == 0

And if you assume Sin[k L]!=0 and divide by that then you get his third equation

In[6]:= Expand[(C[1] Cos[k L] + C[2] Sin[k L])/Cos[k L]] == 0

Out[6]= C[1] + C[2] Tan[k L] == 0

In[7]:= Eliminate[{(L R)/P + C[1] == 0, -(R/P) + k C[2] == 0}, R]

Out[7]= C[1] == -k L C[2] && P != 0

In[8]:= C[1] + C[2] Tan[k L] == 0 /. C[1] -> -k L C[2]

Out[8]= -k L C[2] + C[2] Tan[k L] == 0

In[9]:= Assuming[C[2] != 0, Simplify[%]]

Out[9]= k L == Tan[k L]

Mathematica won' t Solve for or Plot versus expressions, so rename k*L to be kL for a moment

In[10]:= Plot[kL - Tan[kL], {kL, 0, 10}]

Try a few starting x values for FindRoot, all failed to find the "right answer" so do a plot to peek at what we are dealing with.

Out[10]= ...PlotSnipped...

In[11]:= FindRoot[kL == Tan[kL], {kL, 4.5}]

Out[11]= {kL -> 4.49341}

I'd like to reproduce equation 11-42, but I can't see exactly how to do that at the moment.

As you can see, forcing Mathematica to do all this was probably more painful than simply doing the entire thing manually.

Because of the way that Mathematica handles the results, and which was why the author appeared to need to manually force a solution in the form that he wanted, it doesn't appear that it is possible to just let Mathematica determine all the coefficients by giving it all the boundary conditions. That was what we saw earlier. There might be a way to let Mathematica do part of this and still get a result in the form that you are looking for, but at the moment I don't see how.

Is that enough? I have this feeling that this isn't really a satisfactory answer for you and you almost certainly could have done all this by hand in a fraction of the time it took fiddling with trying to get software to do the task.

 

[hushish]

Hi Bill,

Thanks for taking the time to understand the problem. Now you see why, when we add a distributed load q.x into the mix, it complicates things tremendously. I would love to solve the above problem with a q.x added, as well as for the clamped-clamped beam (attached). Then I would be able to solve for all boundary conditions.

If you add a q.x term into the LHS of the equation for the problem you just did above, what particular solution does Mathematica give for the term on the RHS?

Regards,

 

[Bill Simpson]

I'm a little confused here. I expected your last two scanned images to provide the same level of detail for your next two questions. Do I correctly understand that those images just show what we have already done?

I've tried to go back and simplify my previous attempt at a solution.

In[1]:=(v[x]/.DSolve[{v''[x]+k^2 v[x]==R/(e i)(L-x),v[0]==0,v'[0]==0},v[x],x][[1,1]])==0

Out[1]= (k L R - k R x - k L R Cos[k x] + R Sin[k x])/(e i k^3) == 0

Try to automatically separate into a+b*Cos[k x]+c*Sin[k x]

In[2]:= Expand[%]

Out[2]=(L R)/(e i k^2)-(R x)/(e i k^2)-(L R Cos[k x])/(e i k^2)+(R Sin[k x])/(e i k^3)==0

Which after manual substitutions for P=e i k^ 2 and manipulation is

R/P(L-x)-(R L)/P Cos[k x]+R/(P k)Sin[k x]==0

You must justify all these assumptions

In[3]:= Assuming[P != 0 && R != 0 && k != 0, 
 Simplify[R/P(L-x)-(R L)/P Cos[k x]+R/(P k)Sin[k x]==0/.x->L]]

Out[3]= k L Cos[k L] == Sin[k L]

Mathematica doesn't easily Solve for or Plot for expressions so manually replace k*L with kL for this step

In[4]:= Plot[{kL Cos[kL], Sin[kL]}, {kL, 0, 12}]

Out[4]= ...PlotSnipped...

In[5]:= FindRoot[kL Cos[kL] == Sin[kL], {kL, 5}]

Out[5]= {kL -> 4.49341}

If I assume that your next question is this

In[6]:= (v[x] /. DSolve[{v''[x] + (k^2 + q x) v[x] == R/(e i) (L - x), 
      v[0] == 0, v'[0] == 0}, v[x], x][[1, 1]]) == 0

then the result is screens filled with AiryAi, AiryBi, AiryAiPrime, AirBiPrime, BesselJ, HypergeometricPFQ, Gamma and radicals and contains over 8000 individual variables, operators, function names, etc.

If I try this instead

In[7]:= Assuming[q > 0 && k > 0 && R > 0 && L > 0 && x > 0, (Simplify[
    v[x] /. DSolve[{v''[x] + (k^2 + q x) v[x] == R/(e i) (L - x), 
        v[0] == 0, v'[0] == 0}, v[x], x][[1, 1]]]) == 0]

then the result contains over 2000 individual variables, operators, function names, etc.

If I try the corresponding next step and hope for a miracle

In[8]:= Assuming[q > 0 && k > 0 && R > 0 && L > 0 && x > 0, Simplify[% /. x -> L]]

then the result is the same size. I do not understand what information I might be able to provide it to be able to turn this into something that can be understood or used.

If I try this instead

(v[x] /. Assuming[p > 0 && q > 0 && r > 0 && L > 0 && x > 0, 
    Simplify[DSolve[{v''[x]+(p+q x)v[x]==r(L-x),v[0]==0,v'[0]==0}, v[x], x][[1,1]]]])==0

I think I had the result down below 1500 variables, operators, ... but can't seem to reproduce that now.

In that, with the assumption q>0, I see lots of (-q)^ 5/3, (-1)^ 1/3, etc. which I think is Mathematica trying to be careful to not make an error with radicals of negatives. Can you provide information about which root should be chosen in these cases? That might help simplify the problem somewhat.

Have I misunderstood something? If you can provide information to correct my errors or to substantially simplify the problem then there might be a chance. In that case please provide the exact specific equation and assumptions that you want used.

Otherwise I don't think I see any chance of a symbolic solution that can be used. In that case you might be able to construct a numerical solution using your tools, try to verify if that is correct and use that for your problems in the future.

 

[hushish]

Hi Bill,

Again, thanks for the effort.

All your assumptions are correct. The only thing I can add is that Matlab yields the following more-compact result, but I am clueless as to what Int is (it's not integrate):

y =
AiryAi(-(q*E*I)^(1/3)*(P+q*x)/q)*C2+AiryBi(-(q*E*I)^(1/3)*(P+q*x)/q)*C1-R*E*I*pi*(Int(AiryBi(-(q*E*I)^(1/3)*(P+q*x)/q)*(-L+x),x)*AiryAi(-(q*E*I)^(1/3)*(P+q*x)/q)-Int(AiryAi(-(q*E*I)^(1/3)*(P+q*x)/q)*(-L+x),x)*AiryBi(-(q*E*I)^(1/3)*(P+q*x)/q))/(q*E*I)^(1/3)

I have tried various means to arrive at an answer, even substituting the expanded Airy functions, but to no avail. What surprises me is that William F. McCombs seems to have solved this, but infortunately does not give details in his work: "Engineering Column Analysis", 2004. I attach his method.

If I ever find out how to solve it, I will let you know in this thread, and vice versa,

Thanks again,

hushish

 

[AlephZero]

I am trying to recreate the attached graph from Michael Niu, Airframe Stress Analysis and Sizing.

Most likely, that graph (and your other attachment) were made using numerical methods, not by trying to find analytic solutions of a differential equation. This is a trivial problem to model with finite element analysis.

Of course the numerical results will not be "perfectly accurate", but the simplifying assumptions made in this idealized math model of buckling could easily introduce "real world" errors of 50% or more. Getting a mathematically exact solution to the differential equation is not very important compared with that.