Rouben Rostamian

MaplePrimes Activity

These are replies submitted by Rouben Rostamian

After deriving equation (22), the authors identify three special cases for which solutions are available in the literature.  Let's look at their Case 1, that is, a/L=0 and mass_ratio=0.  The equation then should reduce, and they assert that it does, to the classical clamped-free cantilever model.  But I don't see that.

Setting a/L=0 and mass_ratio=0 in equation (22), I get 0=0.  (I did that by hand; no need for Maple for that.)  But that's not correct.  The equation should reduce to cos(βL) cosh(βL) = −1.

Can you verify my calculation?  If I am correct, then there must be a typo in equation (22).  That needs to be fixed before going ahead with the rest of the calculations.


@Neel What you have here is so far removed from what would be the right approach to the solution of the problem, that it leaves little to salvage.  In my earlier message I provided a general outline for solving the forced vibration problem.  What you have in your worksheet does not follow that outline at all.  The entire worksheet needs to be discarded and started over again.

But you should not feel bad about that.  My outline assumes some familiarity with orthogonal basis functions generated by solutions of boundary value problems, and how to expand an arbitrary function in such a basis.  It appears that you have not seen that material, so I wouldn't blame you for that.

I can provide the somewhat long and messy solution to the problem which you are attempting to solve, but without the proper mathematical background it won't be understandable.

Why are you interested in this problem?  If it is a student project, you should really begin with learning the basic mathematical ideas behind what is called eigenfunction expansions. Furthermore, the beam problem is one or two steps more advanced than a beginner problem in this subject, which would be:

First: The heat equation:
w_t= w_xx + f(x,t),
w(0,t)=0,   w(L,t)=0,
w(x,0) = g(x),

Second: The wave equation:
w_tt = w_xx + f(x,t),
w(0,t)=0,   w(L,t)=0,
w(x,0) = g(x),   w_t(x,0) = h(x).

You need to be able to confortably handle these before attempting the beam equation.  Do you have a textbook that shows how to solve these problems?  Learning that will be the first in your progression.  By the way, notice the initial conditions g(x), and h(x).  You will need to specify initial conditions to solve your beam problem as well.  I didn't see initial conditions in your worksheet.

@panke Specify the version of your Maple when you ask a question.  There is a pull-down menu just before you submit your question in which you enter the version.  Please don't ignore it.

The phase space of a 5th order ODE is five-dimensional.  The current version of Maple does three-dimensional pictures quite well.  You will have to wait for a future release of Maple to do plots in five dimensions.

@vercammen If you know about Maple's procs, then a better way would be through a proc, as in

make_prob := proc({n:=10, alpha:=-60,  beta:=20, P_val:=4})
    local Qs := alpha + beta*P;
    printf("There are %a identical firms.  "
           "Supply is given by Qs = %a.  "
           "When P = %a then quantity supplied is equal to %a.  "
           "Calculate the ...\n", n, Qs, P_val, eval(Qs, P=P_val));
end proc:

Then you may call the proc any number of times with the desired arguments.  Missing arguments take on default values as indicated in the first line of the proc's definition.  For example:

make_prob();   # all default args
There are 10 identical firms.  Supply is given by Qs = -60+20*P.  When P = 4 then quantity supplied is equal to 20.  Calculate the ...

make_prob(n=40, P_val=5);   # some optional args
There are 40 identical firms.  Supply is given by Qs = -60+20*P.  When P = 5 then quantity supplied is equal to 40.  Calculate the ...


@Neel There is no general way for that.  Each case needs to be treated on its own, but the steps are pretty much always the same..  So if you see what I did for the clamped-clamped case, you may do the same to analyze the clamped-free case (a cantilever beam).

Here is a pointer to simplify your calculations just a bit.

In what I presented earlier, I went along with general solution of the ODE produced by Maple as a linear combination of exponential and trigonometric functions with coefficients _C1 through _C4.  For our purposes that's not the best form of the general solution.  The first two terms may be replaced with hyperbolic functions, as in _C1*sinh(mu*x/L) + C_2*cosh(mu*x/L).  Since _C1 and _C2 are arbitrary, the changed expression is equivalent to the original.  I suggest that you make that change in the worksheet (manually) and then continue the calculations.  You will see that things work out more smoothly that way.


@Aminaple The idea suggested by Dr. Venkat Subramanian is quite good and I like it better than the solution that I presented earlier.  Here are the details of this alternative approach.


We convert the 2nd order differential equation to a system of 1st order differential
equations by introducing y(t) as the derivative of x(t)

de1 := diff(x(t),t)=y(t);

diff(x(t), t) = y(t)

Then replace "x ' (t) " and x*''*t with y(t) and "y ' (t)" in the differential equation and rearrange:

diff(y(t),t) + 2/25*y(t) + 4*x(t) = (1/25*diff(y(t),t) + 2/3*y(t) + 5/2*x(t)) * x(t-d):
isolate(%, diff(y(t),t)):
de_tmp1 := collect(%, x(t-d));

diff(y(t), t) = (((2/3)*y(t)+(5/2)*x(t))*x(t-d)-(2/25)*y(t)-4*x(t))/(1-(1/25)*x(t-d))

When t < d, the term x(t-d) picks up its values from the history, let call it "h(t)," so the differential

equation really has the form

de_tmp0 := subs(x(t-d)=h(t-d), de_tmp1);

diff(y(t), t) = (((2/3)*y(t)+(5/2)*x(t))*h(t-d)-(2/25)*y(t)-4*x(t))/(1-(1/25)*h(t-d))

Furthermore, when t < 0, we want y(t) to conform to the history, that is,  "y(t)=h ' (t)"  Therefore
"y ' (t) = h ' ' (t)."   We combine the three cases into a single differential equation:

de2 := diff(y(t),t) = piecewise(t<0, (D@@2)(h)(t), t<d, rhs(de_tmp0), rhs(de_tmp1));

de2 := diff(y(t), t) = piecewise(t < 0, ((D@@2)(h))(t), t < d, ((2*y(t)*(1/3)+5*x(t)*(1/2))*h(t-d)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*h(t-d)), ((2*y(t)*(1/3)+5*x(t)*(1/2))*x(t-d)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*x(t-d)))

The initial conditions at t = -d  are determined by the history function:

ic := x(-d) = h(-d), y(-d) = D(h)(-d);

x(-d) = h(-d), y(-d) = (D(h))(-d)

It's time to introduce the history function and the delay:

params := { h = cos,  d = 6/7 };

{d = 6/7, h = cos}

So our system looks like this:

sys := eval({de1, de2, ic}, params);

sys := {diff(x(t), t) = y(t), diff(y(t), t) = piecewise(t < 0, -cos(t), t < 6/7, ((2*y(t)*(1/3)+5*x(t)*(1/2))*cos(t-6/7)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*cos(t-6/7)), ((2*y(t)*(1/3)+5*x(t)*(1/2))*x(t-6/7)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*x(t-6/7))), x(-6/7) = cos(6/7), y(-6/7) = sin(6/7)}

We solve the system and plot the result

dsol := dsolve(sys, numeric):

plots:-odeplot(dsol, [t,x(t)], t=-eval(d, params)..12);

The graph agrees with what I had obtained with the earlier method.






In the forced case you will have
"((&PartialD;)^2w)/((&PartialD;)^( )t^2)+1/(c^(2))((&PartialD;)^(4)w)/((&PartialD;)^( )x^(4))=f(x,t)".

To solve, express both w(x, t) and f(x, t) in terms of the modal functions

"w(x,t)=(&sum;)`alpha__j`(t)*`u__j`(x),  f(x,t)=(&sum;)`beta__j`(t)*`u__j`(x),"

where the coefficients `&beta;__j`(t) are known but `&alpha;__j`(t) are unknown.

Plug this into the PDE, and equate the coefficients of u__j(x).

That will give you an ODE of the form "`alpha__j`' ' + `h__j`*`alpha__j`=`beta__j`" for each j.
Solve those ODEs for `&alpha;__j`(t) and plug the result in the expression

w(x, t) given above.




@Neel What are the boundary condtions?  I cannot tell by reading through your worksheet.

@janhardo Your interpretation of my calculations is quite correct.  Well done!

@AHSAN Even after changing q(y)=j to u(y)=j the problem still makes no sense.  The "bc" and "bc1" in your formulation stand for "boundary condition" but you have not said where the boundaries are.

Go to where you found the problem, read it carefully, then come back with the correct statement.


I looked through your worksheets but was unable to understand them.  It will be very helpful to your readers to have some verbal description of the purpose of the worksheets.  You may want to supply the following information:

  1. It appears that you are solving the PDE 

    with the initial conditions u=(x,0)=sin(Pi*x), diff(u(x,t),t)=0 on the interval (0.1).  Is that correct?
  2. You state the boundary conditions as
    bc1:=u(x,t):        bc2:=u(x,t):
    I don't know that those mean.  Explain.
  3. It appears that you discretize the equation through centered finite differencing in the x direction. Is that correct?
  4. Then you discretize the equations in the time direction in some way.  That does not seem to be consistent with your stated objective of applying the method of lines, which approaches the problem through solving a system of ODEs, but perhaps I have misunderstood.  Explain.
  5. What is the idea behind your non-classical method?
  6. You write "I need to compare the result of the method with the exact solution".  What is the exact solution?
  7. You must know that you may solve your PDE with Maple's pdsolve() (numerically).  You are aware of that, aren't you?
  8. Why have you posted two worksheets?  What is the purpose of each?  Provide a guide.


@janhardo There are neither double integrals nor polar coordinates there.  The area is calculated by simple integration in Cartesian coordinates.

@AHSAN Your bc says q(y) = 0.  Your bc1 says q(y) = j .  Do you want both?


@Fancypants Yes, equation (20) looks good.

1 2 3 4 5 6 7 Last Page 1 of 63