Ali Guzel

Mr. Ali Guzel

80 Reputation

7 Badges

19 years, 111 days

MaplePrimes Activity


These are replies submitted by Ali Guzel

@Ali Guzel  

It is in

The text is from my great Prof Richard H. Enns 2003.

Computer. Algebra Recipes for Classical Mechanics 

Falkland Fiasco in chapter 5.

@acer 

Similar problem in the text that works in Maple 2018
 

restart: with(plots): with(VectorCalculus):  

We assume that g, Omega, V0, and theta are positive, and alpha lies between 0 and (1/2)*Pi.

assume(g>0,Omega>0,V0>0,theta>0,alpha>0,alpha<=Pi/2):

Aliases are introduced.

alias(omega=w,Omega=W,alpha=a):

The angular velocity omega of the Earth at northern latitude alpha  is entered.

w:=<-W*cos(a),0,W*sin(a)>;

Vector(3, {(1) = -W*cos(a), (2) = 0, (3) = W*sin(a)})

(1)

The position vector r of the rocket is given and its velocity v calculated. Here x is to the south, y to the east, and z perpendicular to the surface. These are relative to the noninertial frame attached to the Earth's surface.

r:=<x(t),y(t),z(t)>; v:=diff(r,t);

Vector(3, {(1) = x(t), (2) = y(t), (3) = z(t)})

Vector(3, {(1) = diff(x(t), t), (2) = diff(y(t), t), (3) = diff(z(t), t)})

(2)

The gravitational, Coriolis, centrifugal, and drag forces are entered and the resultant force calculated.

F[gravity]:=<0,0,-g>;

Vector(3, {(1) = 0, (2) = 0, (3) = -g})

(3)

F[Coriolis]:=2*v &x w;

Vector(3, {(1) = 2*(diff(y(t), t))*W*sin(a), (2) = -2*(diff(x(t), t))*W*sin(a)-2*(diff(z(t), t))*W*cos(a), (3) = 2*(diff(y(t), t))*W*cos(a)})

(4)

F[centrifugal]:=w &x (r &x w);

Vector(3, {(1) = -W*sin(a)*(-x(t)*W*sin(a)-z(t)*W*cos(a)), (2) = W^2*cos(a)^2*y(t)+W^2*sin(a)^2*y(t), (3) = -W*cos(a)*(-x(t)*W*sin(a)-z(t)*W*cos(a))})

(5)

F[drag]:=-k*v;

Vector(3, {(1) = -k*(diff(x(t), t)), (2) = -k*(diff(y(t), t)), (3) = -k*(diff(z(t), t))})

(6)

F[resultant]:=F[gravity]+F[Coriolis]+F[centrifugal]+F[drag];

Vector(3, {(1) = 2*(diff(y(t), t))*W*sin(a)-W*sin(a)*(-x(t)*W*sin(a)-z(t)*W*cos(a))-k*(diff(x(t), t)), (2) = -2*(diff(x(t), t))*W*sin(a)-2*(diff(z(t), t))*W*cos(a)+W^2*cos(a)^2*y(t)+W^2*sin(a)^2*y(t)-k*(diff(y(t), t)), (3) = -g+2*(diff(y(t), t))*W*cos(a)-W*cos(a)*(-x(t)*W*sin(a)-z(t)*W*cos(a))-k*(diff(z(t), t))})

(7)

F=ma is applied in the x, y, and z directions.

xeq:=diff(x(t),t,t)=F[resultant][1]; #F=ma in x direction

xeq := diff(x(t), `$`(t, 2)) = 2*(diff(y(t), t))*Omega*sin(alpha)-Omega*sin(alpha)*(-x(t)*Omega*sin(alpha)-z(t)*Omega*cos(alpha))-k*(diff(x(t), t))

(8)

yeq:=diff(y(t),t,t)=F[resultant][2]; #F=ma in y direction

yeq := diff(y(t), `$`(t, 2)) = -2*(diff(x(t), t))*Omega*sin(alpha)-2*(diff(z(t), t))*Omega*cos(alpha)+Omega^2*cos(alpha)^2*y(t)+Omega^2*sin(alpha)^2*y(t)-k*(diff(y(t), t))

(9)

zeq:=diff(z(t),t,t)=F[resultant][3]; #F=ma in z direction

zeq := diff(z(t), `$`(t, 2)) = -g+2*(diff(y(t), t))*Omega*cos(alpha)-Omega*cos(alpha)*(-x(t)*Omega*sin(alpha)-z(t)*Omega*cos(alpha))-k*(diff(z(t), t))

(10)

We use Omega as an expansion parameter and keep only first order terms in Omega in writing x(t)=X0(t)+ Omega X1(t), etc..

x(t):=X0(t)+W*X1(t): y(t):=Y0(t)+W*Y1(t): z(t):=Z0(t)+W*Z1(t):

The above expansions are automatically substituted into xeq, yeq, and zeq.

eq[1]:=expand(xeq); eq[2]:=expand(yeq); eq[3]:=expand(zeq);

eq[1] := diff(X0(t), `$`(t, 2))+Omega*(diff(X1(t), `$`(t, 2))) = 2*Omega*sin(alpha)*(diff(Y0(t), t))+2*Omega^2*sin(alpha)*(diff(Y1(t), t))+Omega^2*sin(alpha)^2*X0(t)+Omega^3*sin(alpha)^2*X1(t)+Omega^2*sin(alpha)*cos(alpha)*Z0(t)+Omega^3*sin(alpha)*cos(alpha)*Z1(t)-k*(diff(X0(t), t))-k*Omega*(diff(X1(t), t))

eq[2] := diff(Y0(t), `$`(t, 2))+Omega*(diff(Y1(t), `$`(t, 2))) = -2*Omega*sin(alpha)*(diff(X0(t), t))-2*Omega^2*sin(alpha)*(diff(X1(t), t))-2*Omega*cos(alpha)*(diff(Z0(t), t))-2*Omega^2*cos(alpha)*(diff(Z1(t), t))+Omega^2*cos(alpha)^2*Y0(t)+Omega^3*cos(alpha)^2*Y1(t)+Omega^2*sin(alpha)^2*Y0(t)+Omega^3*sin(alpha)^2*Y1(t)-k*(diff(Y0(t), t))-k*Omega*(diff(Y1(t), t))

eq[3] := diff(Z0(t), `$`(t, 2))+Omega*(diff(Z1(t), `$`(t, 2))) = -g+2*Omega*cos(alpha)*(diff(Y0(t), t))+2*Omega^2*cos(alpha)*(diff(Y1(t), t))+Omega^2*cos(alpha)*sin(alpha)*X0(t)+Omega^3*cos(alpha)*sin(alpha)*X1(t)+Omega^2*cos(alpha)^2*Z0(t)+Omega^3*cos(alpha)^2*Z1(t)-k*(diff(Z0(t), t))-k*Omega*(diff(Z1(t), t))

(11)

In the following do loop, we produce the equations corresponding to zeroth and first order in Omega .

for i from 1 to 3 do

ode[i]:=collect(lhs(eq[i]),W)=collect(rhs(eq[i]),W): #collect powers of omega

Eqa[i]:=coeff(lhs(ode[i]),W,0)=coeff(rhs(ode[i]),W,0): #equation corresponding to omega^0

Eqb[i]:=coeff(lhs(ode[i]),W,1)=coeff(rhs(ode[i]),W,1); #equation corresponding to omega^1

end do;

ode[1] := diff(X0(t), `$`(t, 2))+Omega*(diff(X1(t), `$`(t, 2))) = (sin(alpha)^2*X1(t)+sin(alpha)*cos(alpha)*Z1(t))*Omega^3+(sin(alpha)^2*X0(t)+sin(alpha)*cos(alpha)*Z0(t)+2*sin(alpha)*(diff(Y1(t), t)))*Omega^2+(2*sin(alpha)*(diff(Y0(t), t))-(diff(X1(t), t))*k)*Omega-k*(diff(X0(t), t))

Eqa[1] := diff(X0(t), `$`(t, 2)) = -k*(diff(X0(t), t))

Eqb[1] := diff(X1(t), `$`(t, 2)) = 2*sin(alpha)*(diff(Y0(t), t))-(diff(X1(t), t))*k

ode[2] := diff(Y0(t), `$`(t, 2))+Omega*(diff(Y1(t), `$`(t, 2))) = (Y1(t)*sin(alpha)^2+Y1(t)*cos(alpha)^2)*Omega^3+(Y0(t)*sin(alpha)^2+Y0(t)*cos(alpha)^2-2*(diff(Z1(t), t))*cos(alpha)-2*sin(alpha)*(diff(X1(t), t)))*Omega^2+(-2*(diff(Z0(t), t))*cos(alpha)-2*sin(alpha)*(diff(X0(t), t))-(diff(Y1(t), t))*k)*Omega-k*(diff(Y0(t), t))

Eqa[2] := diff(Y0(t), `$`(t, 2)) = -k*(diff(Y0(t), t))

Eqb[2] := diff(Y1(t), `$`(t, 2)) = -2*(diff(Z0(t), t))*cos(alpha)-2*sin(alpha)*(diff(X0(t), t))-(diff(Y1(t), t))*k

ode[3] := diff(Z0(t), `$`(t, 2))+Omega*(diff(Z1(t), `$`(t, 2))) = (sin(alpha)*X1(t)*cos(alpha)+cos(alpha)^2*Z1(t))*Omega^3+(sin(alpha)*X0(t)*cos(alpha)+cos(alpha)^2*Z0(t)+2*(diff(Y1(t), t))*cos(alpha))*Omega^2+(-(diff(Z1(t), t))*k+2*(diff(Y0(t), t))*cos(alpha))*Omega-k*(diff(Z0(t), t))-g

Eqa[3] := diff(Z0(t), `$`(t, 2)) = -k*(diff(Z0(t), t))-g

Eqb[3] := diff(Z1(t), `$`(t, 2)) = -(diff(Z1(t), t))*k+2*(diff(Y0(t), t))*cos(alpha)

(12)

The zeroth order equation in the x (south) direction is solved for X0(t), subject to the rocket starting at X0(0)=0 with zero x component of velocity.

sol1:=dsolve({Eqa[1],X0(0)=0,D(X0)(0)=0},X0(t)):

The zeroth order equation in the y (east) direction is solved for Y0(t), subject to the rocket starting at Y0(0)=0 with a velocity component V0*cos(theta).

 

sol2:=dsolve({Eqa[2],Y0(0)=0,D(Y0)(0)=V0*cos(theta)},Y0(t)):

The zeroth order equation in the z  direction is solved for Z0(t), subject to the rocket starting at Z0(0)=0 with a velocity component V0*sin(theta).

sol3:=dsolve({Eqa[3],Z0(0)=0,D(Z0)(0)=V0*sin(theta)},Z0(t)): #solve Eqa[3] for Z0(t)

X0(t), Y0(t), and Z0(t) are now displayed.

X0(t):=rhs(sol1);Y0(t):=rhs(sol2);Z0(t):=rhs(sol3);

X0(t) := 0

Y0(t) := V0*cos(theta)/k-V0*cos(theta)*exp(-k*t)/k

Z0(t) := -exp(-k*t)*(V0*sin(theta)*k+g)/k^2-g*t/k+(V0*sin(theta)*k+g)/k^2

(13)

The above zeroth order solutions are automatically substituted into Eqb[1], Eqb[2], and Eqb[3], which are now shown.

Eqb[1];Eqb[2];Eqb[3];

diff(X1(t), `$`(t, 2)) = 2*sin(alpha)*V0*cos(theta)*exp(-k*t)-(diff(X1(t), t))*k

diff(Y1(t), `$`(t, 2)) = -2*(exp(-k*t)*(V0*sin(theta)*k+g)/k-g/k)*cos(alpha)-(diff(Y1(t), t))*k

diff(Z1(t), `$`(t, 2)) = -(diff(Z1(t), t))*k+2*V0*cos(theta)*exp(-k*t)*cos(alpha)

(14)

Eqb[1] is solved for X1(t), subject to the initial condition that X1(0)=0 and the first order x component of velocity is zero.

sol4:=dsolve({Eqb[1],X1(0)=0,D(X1)(0)=0},X1(t));

sol4 := X1(t) = -(sin(alpha+theta)*V0*k*t+sin(alpha-theta)*V0*k*t+V0*sin(alpha+theta)+V0*sin(alpha-theta))*exp(-k*t)/k^2+V0*(sin(alpha+theta)+sin(alpha-theta))/k^2

(15)

X1(t) is then simplified and the complete x solution to first order in Omega is expressed in XX.

X1(t):=simplify(rhs(sol4));

X1(t) := -V0*(exp(-k*t)*k*t+exp(-k*t)-1)*(sin(alpha+theta)+sin(alpha-theta))/k^2

(16)

XX:=X0(t)+W*X1(t);

XX := -Omega*V0*(exp(-k*t)*k*t+exp(-k*t)-1)*(sin(alpha+theta)+sin(alpha-theta))/k^2

(17)

Eqb[2] is solved for Y1(t), subject to the initial condition that Y1(0)=0 and the first order y component of velocity is zero.

sol5:=dsolve({Eqb[2],Y1(0)=0,D(Y1)(0)=0},Y1(t));

sol5 := Y1(t) = (2*exp(-k*t)*cos(alpha)*g/k+V0*sin(alpha-theta)*(-exp(-k*t)*k*t-exp(-k*t))+2*cos(alpha)*g*t-2*cos(alpha)*g*(-exp(-k*t)*k*t-exp(-k*t))/k-V0*sin(alpha+theta)*(-exp(-k*t)*k*t-exp(-k*t)))/k^2-(V0*sin(alpha+theta)*k-V0*sin(alpha-theta)*k+4*cos(alpha)*g)/k^3

(18)

Y1(t) is then simplified and the complete y solution to first order in Omega is expressed in YY.

Y1(t):=simplify(rhs(sol5));

Y1(t) := (-k*(-1+(k*t+1)*exp(-k*t))*V0*sin(alpha-theta)+(k*V0*(k*t+1)*sin(alpha+theta)+2*cos(alpha)*g*(k*t+2))*exp(-k*t)-V0*sin(alpha+theta)*k+2*cos(alpha)*g*(k*t-2))/k^3

(19)

YY:=Y0(t)+W*Y1(t);

YY := V0*cos(theta)/k-V0*cos(theta)*exp(-k*t)/k+Omega*(-k*(-1+(k*t+1)*exp(-k*t))*V0*sin(alpha-theta)+(k*V0*(k*t+1)*sin(alpha+theta)+2*cos(alpha)*g*(k*t+2))*exp(-k*t)-V0*sin(alpha+theta)*k+2*cos(alpha)*g*(k*t-2))/k^3

(20)

Eqb[3] is solved for Z1(t), subject to the initial condition that Z1(0)=0 and the first order z component of velocity is zero.

sol6:=dsolve({Eqb[3],Z1(0)=0,D(Z1)(0)=0},Z1(t));

sol6 := Z1(t) = -(cos(alpha-theta)*V0*k*t+cos(alpha+theta)*V0*k*t+V0*cos(alpha-theta)+V0*cos(alpha+theta))*exp(-k*t)/k^2+V0*(cos(alpha-theta)+cos(alpha+theta))/k^2

(21)

Z1(t) is then simplified and the complete z solution to first order in Omega is expressed in ZZ.

Z1(t):=simplify(rhs(sol6));

Z1(t) := -V0*(cos(alpha-theta)+cos(alpha+theta))*(exp(-k*t)*k*t+exp(-k*t)-1)/k^2

(22)

ZZ:=Z0(t)+W*Z1(t);

ZZ := -exp(-k*t)*(V0*sin(theta)*k+g)/k^2-g*t/k+(V0*sin(theta)*k+g)/k^2-Omega*V0*(cos(alpha-theta)+cos(alpha+theta))*(exp(-k*t)*k*t+exp(-k*t)-1)/k^2

(23)

The parameter values are entered.

parameters:=g=9.8,V0=500,k=0.01,W=7.27*10^(-5),theta=Pi/3,a=Pi/4:

A function S for subsituting parameters into the various analytic solutions is created and then applied to Y0(t), Z0(t0, XX, YY, and ZZ.

S:=A->evalf(subs(parameters,A)):

Y0(t):=S(Y0(t)); Z0(t):=S(Z0(t)); XX:=S(XX); YY:=S(YY); ZZ:=S(ZZ);

Y0(t) := 25000.00001-25000.00001*exp(-0.1e-1*t)

Z0(t) := -141301.2702*exp(-0.1e-1*t)-0.98e3*t+141301.2702

XX := -2.570333149*exp(-0.1e-1*t)*t-257.0333149*exp(-0.1e-1*t)+257.0333149

YY := 22539.66406-25000.00001*exp(-0.1e-1*t)+94.08072289*(0.1e-1*t+1.)*exp(-0.1e-1*t)+72.70000000*(.1868892204*t+32.54821495)*exp(-0.1e-1*t)+10.07570595*t

ZZ := -141558.3035*exp(-0.1e-1*t)-0.98e3*t+141558.3035-2.570333148*exp(-0.1e-1*t)*t

(24)

We numerically solve for the time tt in seconds for the rocket to hit the ground again. The t range is started at 1 so as to not generate the time that the rocket was fired.

tt:=fsolve(ZZ=0,t=1..500);

tt := 78.40648587

(25)

The southward deflection of the rocket is calculated. For zero air resistance and all other parameters the same, we obtained a deflection of about 100 meters in the text. With air resistance included, the southward deflection is about 1/2 as much.

dx:=eval(XX,{parameters,t=tt});

dx := 47.6775065

(26)

The horizontal range distance traveled  in the y direction is calculated. Without air resistance, we previously obtained a distance of about 22,000 meters. The distance with air resistance is about 2/3 that value.

dy:=eval(YY,{parameters,t=tt});

dy := 13559.29624

(27)

The curve corresponding to the complete solution ([XX, YY, ZZ]) is plotted in gr1 and colored blue. This curve includes gravity, the Coriolus force, the centrifugal force, and air drag.

gr1:=spacecurve([XX,YY,ZZ],t=0..tt,axes=normal,color=blue,thickness=2):

The curve corresponding to Omega = 0  is plotted and colored red.  The Coriolus and centrifugal forces are zero here.

gr2:=spacecurve([X0(t),Y0(t),Z0(t)],t=0..tt,axes=normal,color=red,thickness=2):

The two graphs are displayed together in a 3-dimensional plot which may be rotated by dragging with the mouse. Labels have been added.

display({gr1,gr2},labels=["x (south)","y(east)","z"],tickmarks=[3,3,3],orientation=[-25,60]);

 

 


 

Download 05-S06.mws

@acer 

I have these versions:

Maple 5, Maple 7, Maple 2018, Maple 2021.

The output is from Maple 2026 trial edition

@nm 

Thank you a lot.

@acer 

Thank you my dear prof.

Meanwhile, the Turkish sales representative said that there are no sales of maple versions before 2018

@acer 

I won't respond them.....................

@acer 

lesson22___Suspended_Wire.zip

I agree to you Sir. I would add to your list Prof. Robert Lopez too.

@acer 

I would love to buy you complete lunch, Sir

@acer 

Thank you very much.

Beautifully done. Clean without any errors...

TOP Professor of Mathematics and Maple. Long live Sir

I may have been misunderstood.

The users here, especially Mr. Acer, are very smart and resourceful.

If Prof.Dr. Meade had time, he would have updated the document he produced.

I'm sure of that!

Sincerely,

@Ali Guzel  

I once asked for update of worksheets to Prof Meade. He gave me updated docs except one.  

So, worksheet updates not always easy!!

@Ali Guzel  

I may buy maple 8 since Andre Heck's (2003) book is just excellent..

@acer   

here is Andre Heck's reply to me:

Dear Ali,

On https://staff.fnwi.uva.nl/a.j.p.heck/Maplebook/answers3.html
the link  to the Windows zip-archive works and provides worksheets with all answers to exercises

on https://staff.fnwi.uva.nl/a.j.p.heck/Maplebook/sheets3.html
the link  to the Windows zip-archive works and provides worksheets with all examples in the book

some codes do not run maple 2018.. esp. chp 17

1 2 3 4 Page 1 of 4