Rouben Rostamian

Solution...

Answered: Rouben Rostamian 7711

July 07 2020

1 2

>

restart;

>	with(LinearAlgebra):

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 . 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 :=
a[2]+a[4]=0,
a[1]+a[3]=0,
-b[1]*sin(mu)-b[2]*cos(mu)+b[3]*sinh(mu)+b[4]*cosh(mu)=0,
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 where the matrix and the right-hand side
(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));

(1)

Let us verify that this reduces to the classic in the

case of a cantilever equation which corresponds to .

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

(2)

Nope, that didn't work!

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

order and therefore plugging into it is not the right things to do. Instead,

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

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

-(8/3)*mu^3*(cosh(mu)*cos(mu)+1)

(3)

and there we have the expected equation .

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

in Figure 2 which corresponds to the choice of . 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

>

Download mw.mw

Finding the angle α between P0P1 a...

Answered: Rouben Rostamian 7711

July 06 2020

0 1

Finding the angle α between P₀P₁ and P₀P₂ is a simple matter of applying the law of cosines in the triangle P₀P₁P₂, which is expressed as
||P₂ - P₁||² = ||P₁ - P₀||² + ||P₂ - P₀||² - 2 ||P₁ - P₀|| ||P₂ - P₀|| cos(α),
whence
cos(α) = [ L₁₂² - L₁₀² - L₂₀² ] / [ 2 L₁₀ L₂₀ ].

Is this what you need?...

Answered: Rouben Rostamian 7711

July 05 2020

0 2

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.

restart;
f:=a*u[]*u[1]*u[2]+b*u[1]*u[2]+c*v[]*v[1]*u[2]+d*u[]*v[1]*v[2]+d*u[]*v[1]*u[2]+e*u[]*v[];
    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[]
collect(f,{u[1],u[2],v[1],v[2]},distributed);
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[]]

Maybe this?...

Answered: Rouben Rostamian 7711

July 03 2020

2 3

restart;
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 ...

Needs better care...

Answered: Rouben Rostamian 7711

July 02 2020

1 0

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

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

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

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

Answered: Rouben Rostamian 7711

June 30 2020

2 1

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

Corrections...

Answered: Rouben Rostamian 7711

June 29 2020

3 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};
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);
Then you may plot the solution h1(z) through
plots:-odeplot(dsol, [z, h1(z)], z=0..1);
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.
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!

Needs a refined initmesh...

Answered: Rouben Rostamian 7711

June 28 2020

2 0

>

restart;

>	kernelopts(version);

>	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))-2f_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

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

Download mw.mw

Answer...

Answered: Rouben Rostamian 7711

June 27 2020

3 22

>

restart;

>	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 :

>	eval(pde, w(x,t)=u(x)cos(omegat)); 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 is not convenient. It's better to work with dimensionless quantities.

Therefore we introduce a new unknown, , related to 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

>	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:

>	dsolve(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 , but of course
Maple will pick the trivial solution . The get around that, we temporarily suspend the
boundary condition and instead substitute a boundary condition , where
is an unspecified symbol. Applying the four boundary conditions, we obtain four equations in

the four unknowns 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 , compute the derivative, and set the derivative to zero
in order to assert the boundary condition which we had left out.

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

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 and therefore nonzero.
We conclude that

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

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

Let's verify that the right-hand side of this equation is

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

Okay, therefore eq_tmp is equivalent to:

>	eq := 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=(2n+1)Pi/2), n=1..5);

Call the roots .. For each root there is a mode function given by

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

-(1/2)*a*L^2*(-cos(mu[j])+sin(mu[j])+exp(-mu[j]))*exp(mu[j]*x/L)/(mu[j]^2*(exp(mu[j])-exp(-mu[j])-2*sin(mu[j])))+(1/2)*a*L^2*(exp(mu[j])-cos(mu[j])-sin(mu[j]))*exp(-mu[j]*x/L)/((exp(mu[j])-exp(-mu[j])-2*sin(mu[j]))*mu[j]^2)+(1/2)*a*L^2*(exp(mu[j])-2*cos(mu[j])+exp(-mu[j]))*sin(mu[j]*x/L)/(mu[j]^2*(exp(mu[j])-exp(-mu[j])-2*sin(mu[j])))-(1/2)*a*L^2*cos(mu[j]*x/L)/mu[j]^2

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

multiplicative constant. Also for each there corresponds a frequency 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.

>

Download mw.mw

Solution...

Answered: Rouben Rostamian 7711

June 27 2020

4 1

>

restart;

>	d := 6/7; # the delay

Plot the history function

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

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

>	diff(x(t),t,t) + 2/25diff(x(t),t) + 4x(t) = (1/25diff(x(t),t,t) + 2/3diff(x(t),t) + 5/2x(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 , 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));

Download mw.mw

Solution...

Answered: Rouben Rostamian 7711

June 25 2020

2 4

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

>

restart;

>	C1 := x^2 + y^2 = L^2; C2 := (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;

-(1/8)*L^2*(4*arcsin((1/2)*L/R)*R^2-2*Pi*R^2+L*(-L^2+4*R^2)^(1/2))/R^2

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

(1/8)*(2*R^4*Pi+4*R^4*arcsin((1/2)*(L^2-2*R^2)/R^2)+L^3*(-L^2+4*R^2)^(1/2)-2*L*(-L^2+4*R^2)^(1/2)*R^2)/R^2

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

(1/2)*(-2*arcsin((1/2)*L/R)+Pi)*L^2-(1/2)*L*(-L^2+4*R^2)^(1/2)+(arcsin((1/2)*(L^2-2*R^2)/R^2)+(1/2)*Pi)*R^2

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

Does not exist...

Answered: Rouben Rostamian 7711

June 23 2020

0 1

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

Corrections...

Answered: Rouben Rostamian 7711

June 22 2020

4 3

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:

The mathematical π is written Pi (with uppercase P) in Maple.
The gamma function is written GAMMA (all uppercase) in Maple.
Call your differential equations de1 and de2.
In the line that defines M1, change "=" to ":=".
Enter the initial conditions as
ic := U1(0)=M1, phi(0)=0, D(phi)(0)=0.001;
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.

The longitudinal motion of a rod...

Answered: Rouben Rostamian 7711

June 16 2020

3 1

Here is my paraphrasing of the problem's statement.

The longitudinal vibration of a homogeneous one-dimensional rod of fixed cross-sectional area , modulus of elasticity , and density (mass per unit volume) of , satisfies the PDE where is the displacement of the point at at time , and the subscripts and indicate derivatives. It is customary to let , and write the PDE in the form of the standard wave equation . The strain at any cross-section is , and therefore the stress is , and the force acting on that cross-section is .

Consider such a rod of length , one end of which is fixed at , and the other end is free to move along the axis. The rod's mass is We attach a block of mass to the moving tip of the rod and set it into motion. We wish to model the rod's motion. The initial conditions are 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 of the rod is given by
u(x, t) = r*L*V*(sum(K__j^2*sin(`μ__j`)*sin(`μ__j`*x/L)*sin(c*`μ__j`*t/L)/`μ__j`, j = 1 .. infinity))/c ,
where K__j^2 = 1/(r*sin(`μ__j`)^2+1/2*(1-sin(2*mu[j])/(2*mu[j]))) ,
and the are the solutions of the equation .
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: