> |
with(DEtools):with(plots):
alias(y=y(t), phi=phi(t), y0=y(0),p0=phi(0), yp0=D(y)(0),pp0=D(phi)(0));
eq1 := .10*Diff(y,`$`(t,2))*cos(phi)+.5025*Diff(phi,`$`(t,2))+.490*cos(phi) = PIECEWISE([.2000000000+.559628129599999968*t+.161487481600000010*t^3, t < .5],[.1596281296+.680743740799999997*t+.242231222385861644*(t-.5)^2-1.60743740800000001*(t-.5)^3, t < 1],[.9826030927-.282603092700000002*t-2.16892488954344652*(t-1)^2+3.06826215099999988*(t-1)^3, t < 1.5],[.6254970546-.150331369699999995*t+2.43346833578792321*(t-1.5)^2-2.26561119299999980*(t-1.5)^3, t < 2],[-.5178571430+.583928571299999977*t-.964948453608247214*(t-2)^2+3.99418262299999994*(t-2)^3, t < 2.5],[-5.336542708+2.61461708300000017*t+5.02632547864506574*(t-2.5)^2-10.9111192900000002*(t-2.5)^3, t < 3],[4.027190721-.542396906999999984*t-11.3403534609720165*(t-3)^2+12.8502945499999992*(t-3)^3, t < 3.5],[8.757603092-2.24502945500000006*t+7.93508836524300420*(t-3.5)^2-5.29005890999999994*(t-3.5)^3, otherwise]);
eq2 := 2*Diff(y,`$`(t,2))+.5e-1*cos(phi)*Diff(phi,`$`(t,2))-.5e-1*Diff(phi,t)^2*cos(phi)+5*Diff(y,t)+y+y^3 = PIECEWISE([.3000000000+.401389911599999982*t-.555964654000000014e-2*t^3, t < .5],[.3013899116+.397220176799999990*t-.833946980854194387e-2*(t-.5)^2-.932201767300000039*(t-.5)^3, t < 1],[.8902706185-.310270618499999984*t-1.40664212076583217*(t-1)^2+2.61436671600000015*(t-1)^3, t < 1.5],[.342065540e-1+.243862297300000003*t+2.51490795287187030*(t-1.5)^2-2.40526509599999994*(t-1.5)^3, t < 2],[-1.059642857+.954821428699999974*t-1.09298969072164964*(t-2)^2+1.16669366700000010*(t-2)^3, t < 2.5],[-.642129970+.736851987999999958*t+.657050810014727760*(t-2.5)^2-2.66150957300000002*(t-2.5)^3, t < 3],[3.206688144-.602229381300000033*t-3.33521354933726099*(t-3)^2+5.07934462299999990*(t-3)^3, t < 3.5],[1.347770617-.127934461999999998*t+4.28380338733431554*(t-3.5)^2-2.85586892500000022*(t-3.5)^3, otherwise]);
G:=dsolve({eq1,eq2,y0=0,p0=0,yp0=0,pp0=0.1},[y,phi],'numeric'):
print(" Loi giai so bang phuong phap RUNGE - KUTTA ");
for i from 0 to T do print(G(i)); od;
yy:=t-> rhs(G(t)[2]):
pp:=t-> rhs(G(t)[4]):
yyp:=t->rhs(G(t)[3]):ppp:=t->rhs(G(t)[5]):plot(yy,0..n*T,color=red,thickness=3,title=`tung do y(t)`);
plot(pp,0..n*T,color=blue,thickness=3,title=`goc phi phi(t)`);
plot(yyp,0..n*T,color=green,title=`daohamtungdo y'(t)`);
plot(ppp,0..n*T,color=black,title=`daohamgocphi phi'(t)`); |