JudithRolfe

16 Reputation

2 Badges

16 years, 320 days

MaplePrimes Activity


These are answers submitted by JudithRolfe

My version still doesn't work.  pasted below..

 

restart:
> _Env_dsolve_nowarnstop := true:
>
> #printlevel:=4;
> with(plots):
> with(DEtools):
> traj1:=proc(A,g,cd,spin,phi,print)
> local par1,par2,eq1,eq2,w,ics,pl1,tt,stp,p1,col:
> global sol,a,CL,rho,v,m:
> _Env_dsolve_nowarnstop := true:
> v:=sqrt((diff(x(t),t)^2+(diff(z(t),t)^2))); #sets up the equations
> par1:=(1/2)*A*cd*rho*v;
> par2:=(1/2)*A*CL*rho*v;
> eq1:=m*(diff(x(t),t,t))+par1*diff(x(t),t)+par2*diff(z(t),t)=0;
> eq2:=m*(diff(z(t),t,t))+m*g+par1*diff(z(t),t)-par2*diff(x(t),t)=0;
> #cd:=0.25;
> rho:=1.225;m:=0.05:
> #g:=9.81;
> a:=0.02;
> #A:=Pi*a^2;
> v:=70;
> w:=(Pi*2500/30);
> stp:=[z(t)+10^(-8)];
> #print(A,g,cd,spin);
> #phi:=Pi/4; #sets up the parameters
> ics:=x(0)=0 ,z(0)=0,D(x)(0)=v(0)*cos(phi), D(z)(0)=v(0)*sin(phi); #sets up the initial conditions
> if spin='no' then CL:='0'; #controls the title
> tt:=("without spin");col:=(blue);
> elif A=0 then tt:=("A=0");col:=(black);
> else CL:=0.25+0.1*log(a*w/v);
> tt:=("with spin");col:=(green);
> end if;
> sol:=dsolve({eq1,eq2,ics},{x(t),z(t)},type=numeric,stop_cond=stp);
> if print=yes then
> pl1:=odeplot(sol,[x(t),z(t)],0..12,colour=black,legend=tt,colour=col,title=typeset("Phi= ", phi, "."));
> display(pl1);
> end if;
> end proc:
>
>
To test the procedure - I have asked for 3 separate graphs - with spin, without spin and with A=0;


> traj1(Pi*0.02^2,9.81,0.25,yes,Pi/4,yes);
> traj1(Pi*0.02^2,9.81,0.25,no,Pi/4,yes);
> traj1(0,9.81,0.25,yes,Pi/4,yes);
> traj1(0,9.81,0.25,yes,Pi/4,no);
Warning, cannot evaluate the solution further right of 7.6606639, stop condition #1 violated

any ideas please?

J

Page 1 of 1