@Neel See if the writeup in the attached PDF answers your questions.

By the way, in your line of work you will need to go beyond writing mathematics by hand.  Look into learning LaTeX if you haven't done so already.  My writeup is done in LaTeX.

steady-state-oscillations.pdf

Note added later:
I spotted a few (harmless) typos in the previous PDF.  Here is one with corrections.

steady-state-oscillations2.pdf

You are suffering from culture shock.  You must already know that when in a foreign country, speaking louder in your own language does not make you any better understood.  Learn the local language and customs instead of complaining that they do not agree with yours.

@Neel In equation (9.B.20) you have the force expressed as a delta function.  The picture below that equation also indicates that you do mean a delta function.  But now you are saying q(x,t)= Fo*(1-(x/L))* cos(omega*t) without a delta function.  You need to decide which one you want.

@Neel I looked quickly through your PDF.  The calculations up to (9.B.19) appear to be correct, but after that point I felt lost.

In (9.B.22) we have the equation of motion:

ρ A w_tt + E I w_xxxx = q(x,t)

which is correct, but we have no clear statement of what q(x,t) is.  You need to say
q(x,t) = something.

What is that something?

The boundary conditions
w(0,t) = w_xx(0,t) = w(L,t) = w_xx(L,t) = 0.
for pinned ends are good.

I also don't see what the initial condition are.  You need to specify the initial displacement and initial velocity as

w(x,0) = something
w_t(x,0) = something

Can you tell me what those "something"s are?

@Earl   The Harvard Weekly Problems is a good collection math and physics puzzles.  The one you are referring to is Problem #15.  Beware though that the shape discussed in that problem is the 3D solid obtained by rotating your lamina about the horizontal axis.

@acer  As you have observed, computing  int(int(x,y=-f(x)..f(x)),x=0..1) without the method=meijerg option produces an answer in terms of EllipticE and EllipticK functions with imaginary arguments:

int(int(x, y=-f(x)..f(x)), x=0..1);      
                                     
                      -4/5 EllipticK(I) + 4/5 EllipticE(I)

As we have seen in my original replay, that answer then gets reduced to a real value after applying convert and simplify.

Interestingly, if we change to polar coordinates by letting x=r*cos(t), y=r*sin(t), the equation of the lamina's boundary changes to r = sqrt(cos(t)), where t goes from −π/2 to π/2, and the integration produces an answer in terms of EllipticE and EllipticK functions with real arguments:

int(int(r*r*cos(t), r=0..sqrt(cos(t))), t=-Pi/2..Pi/2);                       
                                   1/2                        1/2
                    1/2           2            1/2           2
              -2/5 2    EllipticK(----) + 4/5 2    EllipticE(----)
                                   2                          2

simplify(convert(%, hypergeom));                                              
                                   1/2           2
                                  2    GAMMA(3/4)
                              2/5 ----------------
                                         1/2
                                       Pi

That is not as good as the direct answer obtained with the method=meijerg option, but at least it has no imaginary parts.

@mmcdara That's an interesting approach.  I wouldn't have thought of applying statistical tools to solve a mechanics problem.

I have converted your "Reply" to an "Answer".  Thumb up!

@Joe Riel That's good. Thanks!

@vv If we could distribute ln over a product, it would convert it to a sum and then the thing will simplify.  But I don't know how to distribute ln even in the simplest case of ln(b/a).  How does one get ln(b) - ln(a) out of this?

@PhD_Wallyson I have made a few corrections to the worksheet and now you can see the modal shapes.  Your interface conditions appear to be correct but I have not examined them carefully.  Be sure to check!  Also be sure that the values of the constants are provided in consistent units.

Download: corrections.mw

@PhD_Wallyson It is possible, in principle, to program Maple to produce a table like that, but I don't think that's worth the effort.  It would be quicker to run your worksheet 10 times with the desired values of the lengths and copy the results into your thesis.

@PhD_Wallyson In this latest worksheet, you handle the symbolic coefficients correctly.  That's good.  Continue doing it that way.

The interface conditions, however, don't look correct.  You have:

(D@@2)(X[1])(L[1]) = -(kr)*(D(X[1])(L[1]))+(D@@2)(X[2])(L[1]),     # equilibrium condition bending moment at x=L[1]

A beam's bending moment is EI X''.  Shoudn't the interface condition have EI in it?

One more thing.  In the worksheet image that you posted today, we see that the μ values are calculated incorrectly.  We even have a negative μ!  I expect that once you correct the missing EI issue, that problem will go away.  And then the plotting of the mode shape will work as before.

Yet another comment: The legnths L1, L2, L3 in the worksheet are in centimeters.  Be sure that kt, kr, E, I, A, etc., are in compatible units.

@PhD_Wallyson Your boundary conditions involve numbers such as 1.422*10^4, 4.881*10^9, etc.  That's makes it impossible for me to verify their correctness.  Please write them symbolically in terms of the parameters as E, I, etc.  Then I will be able to track down the errors, if there are any.

See my previous comments on how to handle numerical values in Maple.

@PhD_Wallyson In my previous reply I noted the need to correct the interface conditions at the supports.  You need to make those corrections before you can compare the results with those of the FEM code.

As to handling the numerical constants, look at our correspondence here.  See how the calculations are carried out without numbers.  Numerical values are introduced through:

span_lengths := L[1]=3.5, L[2]=5.0, L[3]=21.5;
characteristic_equation_numeric := eval(characteristic_equation, {span_lengths});

Do you see how that is done?

@Test007 Specifying the basis is the hard part and it varies greatly from one application to the next.  Do you have something specific in mind?

Once you have the basis, calculating the generalized Fourier series does not require any specialized tools.  You just define the Fourier coefficients through

a[n] = int(f(x)*phi[n](x), x=0..L);   # assuming one dimensional interval (0,L)

and then f(x) is expressed as

f(x) = Sum(a[n]*phi[n](x), n=1..infinity);

That's all.