Rouben Rostamian

MaplePrimes Activity

These are answers submitted by Rouben Rostamian



It is clear from the calculations presented in the paper that the quantities

that matter are mu = beta*L, nu = a/L, rho = m__attached/m__beam.  Their equations would have

been much neater if written in terms of mu, nu, rho.  That's what I will do in what
follows and I suggest that you do it this way in your dissertation.   


Here is the system of linear equations (14) through (21):

sys :=
b[1]*(-cos(mu) + mu*rho*sin(mu)) + b[2]*(sin(mu)+mu*rho*cos(mu))
  + b[3]*(cosh(mu)+mu*rho*sinh(mu)) + b[4]*(sinh(mu)+mu*rho*cosh(mu)) = 0,
a[1]*sin(mu*nu) + a[2]*cos(mu*nu) + a[3]*sinh(mu*nu) + a[4]*cosh(mu*nu) = 0,
b[1]*sin(mu*nu) + b[2]*cos(mu*nu) + b[3]*sinh(mu*nu) + b[4]*cosh(mu*nu) = 0,
a[1]*cos(mu*nu) - a[2]*sin(mu*nu) + a[3]*cosh(mu*nu) + a[4]*sinh(mu*nu)
  - b[1]*cos(mu*nu) + b[2]*sin(mu*nu) - b[3]*cosh(mu*nu) - b[4]*sinh(mu*nu) = 0,
-a[1]*sin(mu*nu) - a[2]*cos(mu*nu) + a[3]*sinh(mu*nu) + a[4]*cosh(mu*nu)
  + b[1]*sin(mu*nu) + b[2]*cos(mu*nu) - b[3]*sinh(mu*nu) - b[4]*cosh(mu*nu) = 0:

We write this in matrix form "M*X=B," where the matrix M and the right-hand side B 
(which is a vector of zeros) are computed through:

M,B :=GenerateMatrix({sys}, [a[1],a[2],a[3],a[4],b[1],b[2],b[3],b[4]]):

To have nonzero solutions to this homogeneous system, we need the determinant

to be zero. That gives us the equation that determines the system's eigenvalues:

eq := simplify(Determinant(M));



Let us verify that this reduces to the classic cos(mu)*cosh(mu)+1 = 0 in the

case of a cantilever equation which corresponds to rho = 0, nu = 0.

eval(eq, {rho=0, nu=0});



Nope, that didn't work!


Looking for the source of trouble, we see that the equation comes with a factor of

order nu^3and therefore plugging nu = 0 into it is not the right things to do.  Instead,

we divide the equation by nu^3, and then take the limit as nu goes to zero:

eval(eq, rho=0):
limit(%/nu^3, nu=0):



and there we have the expected equation cos(mu)*cosh(mu)+1 = 0.


Now with that puzzle out of the way, we turn to producing the plot

in Figure 2 which corresponds to the choice of rho = .2.  Here it is.
You don't need a for-loop for that.

eval(eq, rho=0.2):
plots:-implicitplot(%, nu=0..1, mu=0..18);


You can obtain Figures 3 through 6 in the same way by varying the value of "rho."





Finding the angle α between P0P1 and P0P2 is a simple matter of applying the law of cosines in the triangle P0P1P2, which is expressed as
||P2 - P1||2 = ||P1 - P0||2 + ||P2 - P0||2 - 2 ||P1 - P0|| ||P2 - P0|| cos(α),
cos(α) = [ L122 - L102 - L202 ] / [ 2 L10 L20 ].

You have f both on the right-hand side and on the left-hand side of the definition of f.  I suppose that on the right-hand side you meant to have e instead of f.  I have fixed that.  The list L is the desired result.

    f := a u[] u[1] u[2] + c u[2] v[] v[1] + d u[] u[2] v[1] + d u[] v[1] v[2] + b u[1] u[2] + e u[] v[]
op([1..-1], %);
L := eval([%], [u[1]=1,u[2]=1,v[1]=1,v[2]=1]);
       (a u[] + b) u[1] u[2] + (c v[] + d u[]) u[2] v[1] + d u[] v[1] v[2] + e u[] v[]
       (a u[] + b) u[1] u[2], (c v[] + d u[]) u[2] v[1],  d u[] v[1] v[2], e u[] v[]
       L := [a u[] + b, c v[] + d u[], d u[], e u[] v[]]


n:= 10: alpha:= -60: beta:= 20:
Qs:= alpha + beta*P:
P_val := 4:
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));

Upon execution, this prints:

There are 10 identical firms.  Supply is given by Qs = 20*P-60.  When P = 4 then quantity supplied is equal to 20.  Calculate the ...


@666 jvbasha You can begin by paying closer attention to your code.  In you have

t(x) := sum(p^i*t[i](x), i = 0 .. L);

Do you see that it defines t in terms of ?  Same with the definition of g.  I have not read the code any further.



plot(x^2, x=0..1, labels=[ xi,  `f '`(xi) ] );

  1. Your differential equations are of the first order.  That means that you may specify initial conditions of the values of the unknowns, but not their derivatives. The solution will determine what those derivatives are.  This is what math says, not maple.  Your specification of initial values of derivatives is illegal.  Remove them.  Keep only ics := {h1(0.0001) = 1, h2(0.0001) = 1, h3(0.0001) = 1};
  2. Your system of equations consisting of eq2, eq3, eq4 is a system of ordinary differential equations (ODEs).  Maple's ODE solver is called dsolve.  You are applying pdsolve which is designed for solving parial differential equations (PDEs).  So change the pdsolve line to
    dsol := dsolve(sys1 union ics, numeric);
  3. Then you may plot the solution h1(z) through
    plots:-odeplot(dsol, [z, h1(z)], z=0..1);
  4. When you plot h2 and h3, you will see that they are identical to h1.  That's because your equations and initial conditions are completely symmetrical with respect to h1, h2, h3, and therefore the h1, h2, and h3 are identical.
  5. Additional note: At the top of your worksheet you have:
    with(PDEtools): with(LinearAlgebra): with(DifferentialGeometry): with(VectorCalculus): with(DEtools):
    These are entirely irrelevant to the rest of the worksheet.  Remove!




`Maple 2020.1, X86 64 LINUX, Jun 8 2020, Build ID 1474533`

eq1 := diff(f_n(zeta), zeta$2)-h_n(zeta)*(diff(f_n(zeta), zeta))-f_n(zeta)^2+g_n(zeta)^2 = 0;
eq2 := diff(g_n(zeta), zeta$2)-h_n(zeta)*(diff(g_n(zeta), zeta))-2*f_n(zeta)*g_n(zeta) = 0;
eq3 := 2*f_n(zeta)+diff(h_n(zeta), zeta) = 0;
ics := f_n(0) = 0, f_n(L) = 0, g_n(0) = 1, g_n(L) = 0, h_n(0) = 0;

diff(diff(f_n(zeta), zeta), zeta)-h_n(zeta)*(diff(f_n(zeta), zeta))-f_n(zeta)^2+g_n(zeta)^2 = 0

diff(diff(g_n(zeta), zeta), zeta)-h_n(zeta)*(diff(g_n(zeta), zeta))-2*f_n(zeta)*g_n(zeta) = 0

2*f_n(zeta)+diff(h_n(zeta), zeta) = 0

f_n(0) = 0, f_n(L) = 0, g_n(0) = 1, g_n(L) = 0, h_n(0) = 0

L := 10.0;
dsol := dsolve({eq1,eq2,eq3,ics}, numeric, initmesh = 16):   # default initmesh is 8
   [[zeta,f_n(zeta)], [zeta, g_n(zeta)], [zeta,h_n(zeta)]], color=[red,green,blue]);






pde := diff(w(x,t), x$4) + diff(w(x,t), t$2)/c^2 = 0;

diff(diff(diff(diff(w(x, t), x), x), x), x)+(diff(diff(w(x, t), t), t))/c^2 = 0

We look for solutions of the form w(x, t) = u(x)*cos(omega*t):

eval(pde, w(x,t)=u(x)*cos(omega*t));
de_tmp := expand(%/cos(omega*t));

(diff(diff(diff(diff(u(x), x), x), x), x))*cos(omega*t)-u(x)*omega^2*cos(omega*t)/c^2 = 0

diff(diff(diff(diff(u(x), x), x), x), x)-u(x)*omega^2/c^2 = 0

The coefficient omega^2/c^2 is not convenient.  It's better to work with dimensionless quantities.

Therefore we introduce a new unknown, mu, related to omega through

mu_def := omega = c/L^2*mu^2;

omega = c*mu^2/L^2

Then express the differential equation de_tmp in terms of mu

de := eval(de_tmp, mu_def);

diff(diff(diff(diff(u(x), x), x), x), x)-u(x)*mu^4/L^4 = 0

Here is the general solution of the de:

U := unapply(rhs(%), x);

proc (x) options operator, arrow; _C1*exp(mu*x/L)+_C2*exp(-mu*x/L)+_C3*sin(mu*x/L)+_C4*cos(mu*x/L) end proc

We want to apply the boundary conditions U(0) = 0, U*('0 = 0, U(L) = 0, U')(L) = 0, but of course
Maple will pick the trivial solution U(x) = 0.  The get around that, we temporarily suspend the
boundary condition "U '(L)=0," and instead substitute a boundary condition (diff(U(x), x, x))*0 = a, where a
is an unspecified symbol.  Applying the four boundary conditions, we obtain four equations in

the four unknowns _C1, _C2, _C3, _C4 which we solve:

U(0) = 0, D(U)(0)=0, (D@@2)(U)(0) = a, U(L) = 0:
C_vals := solve({%}, {_C1, _C2, _C3, _C4});

{_C1 = -(1/2)*a*L^2*(-cos(mu)+sin(mu)+exp(-mu))/(mu^2*(exp(mu)-exp(-mu)-2*sin(mu))), _C2 = (1/2)*a*L^2*(exp(mu)-cos(mu)-sin(mu))/((exp(mu)-exp(-mu)-2*sin(mu))*mu^2), _C3 = (1/2)*a*L^2*(exp(mu)-2*cos(mu)+exp(-mu))/(mu^2*(exp(mu)-exp(-mu)-2*sin(mu))), _C4 = -(1/2)*a*L^2/mu^2}

Plug these into the expression for U(x), compute the derivative, and set the derivative to zero
in order to assert the boundary condition "U '(L)=0" which we had left out.

eval(U(x), C_vals):
diff(%, x):
eval(%, x=L) = 0:

a*L*(exp(2*mu)*cos(mu)-2*exp(mu)+cos(mu))/(mu*(exp(2*mu)-2*sin(mu)*exp(mu)-1)) = 0

The denominator of that equation is a multiple of sinh(mu)-sin(mu) and therefore nonzero.
We conclude that

exp(2*mu)*cos(mu) - 2*exp(mu) + cos(mu) = 0;
eq_tmp := isolate(%, cos(mu));

exp(2*mu)*cos(mu)-2*exp(mu)+cos(mu) = 0

cos(mu) = 2*exp(mu)/(exp(2*mu)+1)

Let's verify that the right-hand side of this equation is 1/cosh(mu)NULL

simplify(1/cosh(mu) - rhs(eq_tmp));


Okay, therefore eq_tmp is equivalent to:

eq := cos(mu) = 1/cosh(mu);

cos(mu) = 1/cosh(mu)

This equation has infinitely many roots, as can be seen in this graph:

plot([cos(mu),1/cosh(mu)], mu=0..5*Pi, color=[blue,red], thickness=3, size=[700,300]);

Here are the first few roots:

seq(fsolve(cos(mu)=1/cosh(mu), mu=(2*n+1)*Pi/2), n=1..5);

4.730040745, 7.853204624, 10.99560784, 14.13716549, 17.27875966

Call the roots `μ__j`, j = 1, 2, () .. ()..  For each root there is a mode function given by

subs(C_vals, mu=mu[j], U(x));


where the factors a may be dropped since a mode function is determined modulo a

multiplicative constant.  Also for each `μ__j` there corresponds a frequency `ω__j`  given through

the equation mu_def, that is, "`omega__j`=(c*`mu__j`^(2))/(L^(2))."  You should be all set now for eigenfunction expansion

of the solutions of your PDE.





d := 6/7;  # the delay


Plot the history function

p[0] := plot(cos(t), t=-d..0):

The solution on the interval 0, d picks up the history function (colored green)

diff(x(t),t,t) + 2/25*diff(x(t),t) + 4*x(t)
  = (1/25*diff(x(t),t,t) + 2/3*diff(x(t),t) + 5/2*x(t))*cos(t-d),
  x(0)=1, D(x)(0)=0:
dsol := dsolve({%}, numeric, output=operator, range=0..d):
myx[1] := eval(x, dsol):
mydx[1] := eval(D(x), dsol):
p[1] := plot(myx[1](t), t=0..d):

On each subsequent interval [(i-1)*d, i*d], the differential equation picks up the solution

calculated in the preceding interval (colored green)

for i from 2 to 14 do
sys := diff(x(t),t,t) + 2/25*diff(x(t),t) + 4*x(t)
  = (1/25*diff(x(t),t,t) + 2/3*diff(x(t),t) + 5/2*x(t))*myx[i-1](t-d),
  x((i-1)*d)=myx[i-1]((i-1)*d), D(x)((i-1)*d)=mydx[i-1]((i-1)*d):
dsol := dsolve({sys}, numeric, output=operator, range=(i-1)*d..i*d):
myx[i] := eval(x, dsol):
mydx[i] := eval(D(x), dsol):
p[i] := plot(myx[i](t), t=(i-1)*d..i*d);
end do:
plots:-display(seq(p[j], j=0..i-1));



Here is the solution with no explanations.  See if you can figure it out.


C1 := x^2 + y^2 = L^2;
C2 := (x-R)^2 + y^2 = R^2;

x^2+y^2 = L^2

(x-R)^2+y^2 = R^2

solve(C1, y):
c1 := %[1];


solve(C2, y):
c2 := %[1];


X := solve(c1=c2, x);


A1 := int(c1, x=X..L) assuming L>0, 2*R>L;


A2 := int(c2, x=0..X) assuming L>0, 2*R>L;


A := simplify(2*(A1+A2));


A = Pi*R^2/2:
eval(%, R=1):
fsolve(%, L=0..2);



Such a rational function does not exist!  Is that a trick question?

Since you are not in the habit of posting a worksheet, I won't post a  worksheet either.  I will give my comments as text.

Make the following corrections in your code:

  1. The mathematical π is written Pi (with uppercase P) in Maple.
  2. The gamma function is written GAMMA (all uppercase) in Maple.
  3. Call your differential equations de1 and de2.
  4. In the line that defines M1, change "=" to ":=".
  5. Enter the initial conditions as
    ic := U1(0)=M1, phi(0)=0, D(phi)(0)=0.001;
  6. Set T:=0.1;  beta:=0.6;  k:=3.5;

After these changes you should be able to call Maple's dsolve to solve your system numerically:

dsol := dsolve({de1,de2,ic}, numeric);

To plot U1, phi, and the derivative of phi over the time interval 0 to 30, do

plots:-odeplot(dsol, [x,U1(x)], x=0..30, view=0..1);
plots:-odeplot(dsol, [x,phi(x)], x=0..30);
plots:-odeplot(dsol, [x,diff(phi(x),x)], x=0..30);

You will see that Maple issues a warning, saying that it cannot go beyond 23.53.  If you look closely for the reason for that, you will see that at t=23.53 the value of phi is zero and the value of its derivative is negative, therefore at the next timestep the value of phi is going to turn negative.  But that is problematic since your differential equation contains the 3/2 power of phi.  What is the 3/2 power of a negative number?  Is there a reason why you have that term in your equation?  Perhaps it shouldn't be there.  Or perhaps it should be the absolute value of phi raised to the 3/2 power.  Go back to the origin of the equation and investigate.


Here is my paraphrasing of the problem's statement.

The longitudinal vibration of a homogeneous one-dimensional rod of fixed cross-sectional area A, modulus of elasticity E, and density (mass per unit volume) of rho, satisfies the PDE   rho*u__tt = E*u__xx"," where u = u(x, t) is the displacement of the point at x at time t, and the subscripts x and t indicate derivatives.  It is customary to let c^2 = E/rho, and write the PDE in the form of the standard wave equation u__tt = c^2*u__xx. The strain at any cross-section is u__x, and therefore the stress is E*u__x, and the force acting on that cross-section is E*A*u__x.


Consider such a rod of length L, one end of which is fixed at x = 0, and the other end is free to move along the x axis.  The rod's mass is "M=rho*A*L."  We attach a block of mass m to the moving tip of the rod and set it into motion.  We wish to model the rod's motion. The initial conditions are "u(x,0)=0,"and
"((&PartialD; u)/(&PartialD; t)) ? ()|() ? (t=0)={[[0,        0<=x<L],[V,x=L]]"

This is a somewhat tricky problem since the initial velocity is discontinuous at the boundary. In the attached worksheet I show that the displacement u(x, t) of the rod is given by
u(x, t) = r*L*V*(sum(K__j^2*sin(`&mu;__j`)*sin(`&mu;__j`*x/L)*sin(c*`&mu;__j`*t/L)/`&mu;__j`, j = 1 .. infinity))/c,
where "r=m/(M),"   K__j^2 = 1/(r*sin(`&mu;__j`)^2+1/2*(1-sin(2*mu[j])/(2*mu[j]))),
and the `&mu;__j` are the solutions of the equation mu*tan(mu) = 1/r.
Here is a sample of what the solution looks like (u(x,t) versus x, at various times t):

and here is an animation of the rod itself:


The short vertical lines are equally-spaced markers etched on the rod to help visualize the rod's motion.


Download worksheet:




T := v -> sum(v(i)*x^(i-1),i=1..n);

proc (v) options operator, arrow; sum(v(i)*x^(i-1), i = 1 .. n) end proc

We want to show that T(x)+T(y) = T(x+y), and T(c*x) = c*T(x) for all
n-tuples x and y, and all scalars c. The first one is easy:

combine(T(x)+T(y)) = T(x+y);

sum(x(i)*x^(i-1)+y(i)*x^(i-1), i = 1 .. n) = sum((y(i)+x(i))*x^(i-1), i = 1 .. n)


The second one takes a bit more of effort since we need to tell Maple that
x is an n-tuple while c is a scalar.  We have


sum(c(i)*x(i)*x^(i-1), i = 1 .. n)

which is not what we want since c(i) should be just c.  To fix that, we define c
to be the constant function that assigns the constant c to any t, that is, c = (proc (t) options operator, arrow; c end proc)

eval(T(c*x), c=(t->c));

sum(c*x(i)*x^(i-1), i = 1 .. n)

Now that's better.  Let's verify that T(c*x) = c*T(x):

eval(T(c*x), c=(t->c)) = c*T(x);

sum(c*x(i)*x^(i-1), i = 1 .. n) = c*(sum(x(i)*x^(i-1), i = 1 .. n))




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