@Carl Love Hi Carl, thanks for pointing that out!

I'm not sure yet whether the value makes sense or not, I need to test it with other values first.

However, the shape of the function should not be a half-wave. Based on the experiments, the shape of the function should have around 3-5 half-waves, not just one. Does it mean that I might need the next eigenvalue to show the shape of this function? And is there a way to do that?

Thanks again!

@Carl Love @Rouben Rostamian Hi Carl and Rouben, I tried the method you suggested, but however, the eigenvalue I got is the trivial eigenvalue which is 0 (sigma=0). Is there a way to find the next lowest eigenvalue using this method?

Please find attached the worksheet Longitudinal_ODE_2.mw

Thanks!

@Carl Love Hi Carl, thanks again for your explanation! I appreciate your time and effort to explain it!

One of the reasons I want to get rid of the imaginary part is that the analysis takes too long when the imaginary part occurs (around 5-10 minutes) (but only to get 10^-13*I or 0*I). Some analysis I have done previously, where only the real part exists, it takes seconds to finish it. 

As you mentioned there is an easy way to remove the imaginary part, Could you please elaborate on how to do it cause it might speed up my code.

I'm a bit rusty in my calculus since it's been at least 3 years since I took my calculus class in uni. I didn't expect I might go through ODE in my thesis. Thanks again for the reminder and your explanations!

About the Numeric BVP Solver, as you mentioned that Rouben's method would not work, could you explain more about how I can get my eigenvalues again?

when I do sol_L:=dsolve({deq||(2..6)}union{eq||(2..21)},numeric); (where deq2..6 are my ODE and eq2..21 are my BC and end points interval), I get trivial solutions as expected (Vup1,Vp1,Vp2,Vlp1,Vlp2 = 0). What term should I add to be able to find my eigenvalue? Do I do this? dsolve({deq||(2..6)}union{eq||(1..21)},numeric); (where eq1 is BC involving sigma) and leave sigma as a parameter?

Thanks again, Carl!

@Rouben Rostamian  Hi Rouben, Thanks for your suggestions! I'll try and let you know soon!

@Rouben Rostamian  Hi Rouben, I actually modelled a folded plate, which plates are connected to each other in y direction through continuity equations.

I only apply the buckling (geometric non-linearity) in x direction, while in y direction the plates are assumed to be linear elastic, simplifying the problem similar to a beam.

The y linear equation can be found in "transverse behaviour" in my original spreadsheet. The idea is that I substitute all the behaviour of transverse direction in term of x from the continuity equations, and subs it into longitudinal equations.

The boundary condition is simply supported which means at 0 and L, V of each panels and the bending moment (EI*diff(V,x$2)) are 0.

I did the buckling of a simple beam and it works fine.

Since if I do sol_"L:=evalf(value(dsolve({deq||(2..6)},{Vup1(x),Vp1(x),Vp2(x),Vlp1(x),Vlp2(x)})));" I got solutions (even it's in complex number), is it safe to assume that the problem is in the boundary condition?

Or is it the other way round? (The problem is in the ODE since it yields complex numbers)

Thanks again!

@Rouben Rostamian @mmcdara Also, If I put floating point in the first place as a parameter (0.1 instead of 1/10) it gives me complex number again. But again the imaginary term of the complex number is close to 0. (i.e 0*I, 10^-13*I).

@Carl Love Thanks for the explanation!

@Rouben Rostamian  @mmcdara Hi all

The physical meaning of my model is to find a buckling stress (sigma) (geometric non-linearity) of plates. Basically, what I'm doing in model 1 is that I didn't include the buckling stress (sigma) (geometric non-linearity) in the model, and I found the deflection (Vup1(x),Vp1(x),Vp2(x),Vlp1(x),Vlp2(x),etc) is approximately the same as experiments and simple model. To find the buckling stress (sigma), I need to add the geometric non-linearity which is the last term in my ODE equations (sigma*A*diff(V(x),x$2)) for each panel. After I put the boundary condition, as @mmcdara noticed that the trivial solution for this problem is all V equals to 0.

I need to find what sigma value where V of each panel does not equal to 0 (non-trivial solution).

My professor told me to solve the problem is to create a matrix from the boundary condition, (in the form like 1*c1+3*c2+6*c3=5*sigma*A) and make it a matrix (say matrix B), hence this will become an eigenvalue problem, where I have to find what sigma value to make determinant of matrix B = 0. Then it will be what my buckling stress it.

But somehow I'm getting the complex number when solving the general solutions after I add the geometric nonlinearity (sigma*A*diff(V(x),x$2)). I tried to rewrite the equations in a separate worksheet and solve it simultaneously with the boundary condition as you guys suggest, I'm getting empty set for my sol_L. But if I put sigma:=0, (basically back to model 1), I'm getting results again.

Please find the worksheet attached Longitudinal_ODE.mw

Thanks a lot everyone, I really appreciate your willingness to help!

@mmcdara Hi, as mentioned in my reply before, I need the coefficient of each integration coefficient (after I put the boundary condition) to make a matrix.

I cannot do Vp1(0), because it just gives me Vp1(0) back, (when I type Vp1(0); it shows Vp1(0) as an output) that's why I do subs(x=0,Vp1(x));

Also, as I am a beginner in maple, do you mind explaining what D@@2 does?

Or is there a page that I can read about solving ODE and its boundary condition at the same time?

Thanks a lot!

@Rouben Rostamian  Hi Rouben, the reason why I solve for the general solution first is that I need the coefficient of each integration coefficient (after I put the boundary condition) to make a matrix and solve for an eigenvalue for later use. Currently, I'm just trying to see whether the code is working or not.

For example, after I put the boundary condition and I yield something like eq1:=1*c1+3*c2+6*c3, I need to make a matrix with the first row of [1 3 6].

In fact, after I put the boundary condition, the final form of each function is real number (with 0*I or 10^-12*I) like Vp1(x):= 5*cos(x)+0*I.

Also I'm not sure about how to solve the ODE with boundary condition at the same time actually, do you mind explaining to me how? Or is there a webpage that you can show me?

Thanks again!

@Carl Love Hi Carl, yeah, it was my mistakes. When I saw /tan(Pi/2) my initial reaction is to add epsilon to remove the problem. I should've thought to change it into *cot(Pi/2). 

Thanks again! It definitely taught me a lesson!

@tomleslie Hi Tom, sorry for the late reply. The program is working perfectly fine now. Thanks for your help! Appreciate that :)

@Carl Love Thanks for your feedback, Carl

It will really help me in my future as well.

I appreciate that.

@Carl Love Hi Carl, I've already done what you suggested (remove epsilon and change /tan(Pi/2) into *cot(Pi/2).

However, the erratic behaviour occurs again every time I rerun with the same parameter.

Please find attached the updated worksheet and the screenshotsSP1JH7.2.mw 

@Carl Love Hi Carl,

Thanks for your feedback.

By what I meant the results I want is that the result is similar when comparing to the experimental results.

In that case, I'll redo my model based on your comments too.

Thanks!