# M145L.W99, Week 4, Lab #3
# by ralph abraham, 24 january 1999
# file: ueda, this is a MapleVr5 worksheet 
# 
# we follow Ueda p. 208
# this also plots a Poincare trajectory
# This is the forced Duffing in Ueda form
# Experiment with changes in the parameters (k > 0, B >= 0)
# Define the functions
> k := 0.1:
> B := 12.0:
> xdotexpr := y:
> ydotexpr := -k*y-x*x*x+B*cos(t):
# Define the odes
> myvector1 := subs({x=X(t), y=Y(t)}, xdotexpr):
> myvector2 := subs({x=X(t), y=Y(t)}, ydotexpr):
> eq1 := diff(X(t),t)=myvector1:
> eq2 := diff(Y(t),t)=myvector2:
# 
# Define Poincare first return function
> P := proc() local ini, sol, x1, y1;
> ini := X(0) = args[1][1], Y(0) = args[1][2];
> sol := dsolve( {eq1, eq2, ini}, {X(t),Y(t)}, type=numeric );
> x1 := subs(sol(2*Pi)[2], X(t));
> y1 := subs(sol(2*Pi)[3], Y(t));
> [x1,y1];
> end:
# 
# 
# Iterate the Poincare map
# Trajectory
# 
> N := 10;

                               N := 10

> x0 := 1.0;

                              x0 := 1.0

> y0 := 0.0;

                               y0 := 0

> initpt := [x0,y0];

                          initpt := [1.0, 0]

> traj := array(0..N);

                      traj := array(0 .. 10, [])

> traj[0] := initpt;

                         traj[0] := [1.0, 0]

> for n from 1 to N do
>   traj[n] := P(traj[n-1]);
> od;

          traj[1] := [1.159358846033247, -1.764986337454787]


          traj[2] := [.8159922049208893, -1.968860662312962]


          traj[3] := [.7010815979645595, -1.355111330724975]


          traj[4] := [.8395003431147330, -.9318907193308145]


          traj[5] := [.9706406934679841, -1.165089606694018]


          traj[6] := [.9282161585902112, -1.470289168652357]


          traj[7] := [.8540018184172132, -1.390419580060766]


          traj[8] := [.8621257215559371, -1.223176959267183]


          traj[9] := [.9043173684036405, -1.221102360421068]


         traj[10] := [.9098203905979048, -1.314790851007412]

> with(plots):
> pointplot(traj);

# This is the end of file: ueda03