Question: Maple 9.5 example, in later version

05-2-2.mws

Can you help me with this code?

restart: with(VectorCalculus):

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

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

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

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

(1)

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)

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

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

(3)

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

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

(4)

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

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

(5)

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

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

(6)

eq:=(u,i)->simplify(diff(u(t),t,t)=F[resultant][i]):

xeq:=eq(x,1); yeq:=eq(y,2); zeq:=eq(z,3);

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

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

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

(7)

ic:=x(0)=0,y(0)=0,z(0)=0,D(x)(0)=0,D(y)(0)=V0*cos(theta),D(z)(0)=V0*sin(theta);

ic := x(0) = 0, y(0) = 0, z(0) = 0, (D(x))(0) = 0, (D(y))(0) = V0*cos(theta), (D(z))(0) = V0*sin(theta)

(8)

sol:=dsolve({xeq,yeq,zeq,ic},{x(t),y(t),z(t)},method=laplace):

assign(sol):

f:=u->simplify(expand(u(t))): X:=f(x); Y:=f(y); Z:=f(z);

X := -(1/4)*(Omega^4*V0*sin(theta)*cos(alpha)*(sum(exp(_alpha1*t)/((Omega^2+_alpha1^2)*_alpha1), _alpha1 = RootOf(Omega^2+_Z^2)))+cos(alpha)*g*(sum(exp(_alpha1*t)*_alpha1^2/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))-sin(theta)*cos(alpha)*V0*(sum(exp(_alpha1*t)*_alpha1/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))*Omega^2+(-2*Omega^3*cos(theta)*V0+3*cos(alpha)*Omega^2*g)*(sum(exp(_alpha1*t)/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))+4*(Omega^2*t*V0*sin(theta)-(1/2)*Omega^2*t^2*g-g)*cos(alpha))*sin(alpha)/Omega^2

Y := (1/4)*(-V0*cos(theta)*Omega^2+2*cos(alpha)*Omega*g)*(sum(exp(_alpha1*t)/((Omega^2+_alpha1^2)*_alpha1), _alpha1 = RootOf(Omega^2+_Z^2)))-(1/2)*sin(theta)*cos(alpha)*V0*Omega*(sum(exp(_alpha1*t)/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))+(1/4)*V0*cos(theta)*(sum(exp(_alpha1*t)*_alpha1/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))

Z := (1/4)*(-Omega^4*cos(alpha)^2*V0*sin(theta)*(sum(exp(_alpha1*t)/((Omega^2+_alpha1^2)*_alpha1), _alpha1 = RootOf(Omega^2+_Z^2)))-cos(alpha)^2*g*(sum(exp(_alpha1*t)*_alpha1^2/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))+cos(alpha)^2*sin(theta)*V0*(sum(exp(_alpha1*t)*_alpha1/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))*Omega^2+(2*Omega^3*cos(alpha)*cos(theta)*V0-3*Omega^2*g*cos(alpha)^2)*(sum(exp(_alpha1*t)/(Omega^2+_alpha1^2), _alpha1 = RootOf(Omega^2+_Z^2)))+(-4*Omega^2*t*V0*sin(theta)+2*Omega^2*t^2*g+4*g)*cos(alpha)^2+4*Omega^2*t*V0*sin(theta)-2*Omega^2*t^2*g)/Omega^2

(9)

P:=(u,n)->convert(taylor(u,W=0,n),polynom):

Xexp:=P(X,4); Yexp:=P(Y,4); Zexp:=P(Z,4);  

Error, (in series/sum) unable to compute series

Error, (in series/sum) unable to compute series

Error, (in series/sum) unable to compute series

 

tt:=solve(Zexp=0,t);

tt :=

(10)

 T1:=P(tt[2],1); d[x]:=eval(Xexp,t=T1);  

Error, invalid subscript selector

d[x] := Xexp

 

T2:=P(tt[2],2); d[y]:=P(eval(Yexp,t=T2),2);

Error, invalid subscript selector

d[y] := Yexp

 

d[y]:=collect(d[y],[cos(a),1/g^2,V0^3,W]);

d[y] := Yexp

(11)

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

d[x]:=eval(d[x],evalf(parameters));

d[x] := Xexp

(12)

d[y]:=eval(d[y],evalf(parameters));

d[y] := Yexp

(13)
 

 

Download 05-2-2.mws

Please Wait...