Numerical solution of damped, driven pendulum. Author: David Craig <http://web.lemoyne.edu/~craigda/> Version: 1.9 [5-10-06] To do: (i) Include bifurcation diagram? [The change numsol -> numsol(g) necessary to generate it slows everything else down a bit, so perhaps better in separate file] (ii) Improve procedure for resetting angle x mod 2Pi \342\200\223 in particular, eliminate "connectives" at jumps [plot as multiple plots over [-Pi+delta,+Pi-delta], say?] (iii) Improve titles/legends (iv) More elegant way to define difference of neighboring solutions Dnumsol restart; with(plots):Digits := 15; # numerical precision [default value is 10]; larger is slowerDamped driven oscillator parameters: # default values as in Taylor, chapter 12 w0 := (3/2)*w; # natural frequency; b := w0/4; # damping parameter [Taylor's beta] g := 0.6; # driving parameter [Taylor's gamma] w := 2*Pi; # driving frequency T := 2*Pi/w; # driving period x - angle; y - angular velocity; z - wt Equations of motion in first order form: eq1 := diff( x(t),t ) = y(t); eq2 := diff( y(t),t ) = -2*b*y(t) - w0^2*sin(x(t)) + g*w0^2*cos(z(t)); eq3 := diff( z(t),t ) = w; Initial conditions x0 := -Pi/2; y0 := 0; z0 := 0; ic := x(0)=x0,y(0)=y0,z(0)=z0; list of parameter values for title/legend plist := cat("w=",convert(w,string)," ", "T=",convert(T,string)," ", "w0=",convert(w0,string)," ", "b=",convert(b,string)," ", "g=",convert(g,string), "\n", "x0=",convert(x0,string)," ", "y0=",convert(y0,string) ) ; # options for appearance of plots plotopts := 'axes=BOXED,title=plist,titlefont=[HELVETICA,16]' ; Numerical solution with initial conditions ic [consider method=gear at some point] (The option maxfun=0 eliminates the errors that throw when you try to evaluate to large t [> 90, roughly]. See Help for dsolve,maxfun.) numsol := dsolve( [eq1,eq2,eq3,ic],[x(t),y(t),z(t)],maxfun=0,numeric );Output of dsolve is a procedure that numerically evaluates the functions [x(t),y(t),z(t)] at the time specified by its argument numsol(1); # for example, this is a list of the values of those functions at t=1 numsol(1)[3]; # this is how you access the list results individually - note each list entry is an EQUATION # (see the help for "dsolve,numeric" for details and instructions on options to change this) # this is how you get the actual numerical value of each result - rhs() extracts the right side of an equation rhs( numsol(1)[3] ) ; Resets x to remain in the range -Pi..Pi [ i.e. x mod 2Pi] # subtracts/adds the appropriate multiple of 2Pi when x increases past Pi/decreases past -Pi mod2Pi := x -> x - 2*Pi*floor( (x+Pi)/(2*Pi) ) :Angle vs. time odeplot( numsol, [t,x(t)],0..30,plotopts,numpoints=400); # same plot but with x mod 2Pi # odeplot( numsol, [t,mod2Pi(x(t))],25..30,plotopts,numpoints=400); # plots x mod2Pi Animate trajectory ; only execute when you need it as this is resource intensive odeplot( numsol, [t,x(t)],0..30,plotopts,frames=400);Phase-space trajectory odeplot( numsol, [x(t),y(t)],0..30,plotopts,numpoints=3000); # same plot but with x mod 2Pi # odeplot( numsol, [mod2Pi(x(t)),y(t)],30..60,plotopts,numpoints=1400); Animate phase trajectory ; only execute when you need it as this is QUITE resource intensive - and sometimes twitchy odeplot( numsol, [x(t),y(t)],0..30,plotopts,frames=400);phase-space trajectory in time [3D] odeplot( numsol, [x(t),y(t),z(t)],0..10,plotopts,numpoints=400);Poincare section [NB: Maple also has its own built-in tools for Poincare sections!] N := 60; # number of periods to plot ni := 30; # starting period of section nf := ni + N-1; # ending period of section ti := ni*T ; # starting time of section tf := nf*T; # ending time of section wx := 25; # square plot window bounds \342\200\223 x wv := 25; # square plot window bounds \342\200\223 v pwsq := view=[-v..v,-v..v] ; # square plot window pwrct := view=[-wx..wx,-wv..wv] ; # rectangular plot window pwot := view=[-10..0,15..25] ; # other plot window # add (ni,N) to parameter list in plot title plist := cat(plist,"\n", "ni=",convert(ni,string)," ", "N=",convert(N,string) ) ; # numpoints = N ensures a point is plotted just once per driving period odeplot( numsol, [x(t),y(t)],ti..tf,plotopts,style=point,symbol=circle,symbolsize=5,numpoints=N,pwrct); # same plot but with x mod 2Pi # odeplot( numsol, [mod2Pi(x(t)),y(t)],ti..tf,plotopts,style=point,symbol=circle,symbolsize=5,numpoints=N,pwrct); Animate Poincare section ; only execute when you need it as this is resource intensive \342\200\223 and twitchy? # animation shows additional points at early times I don't yet understand odeplot( numsol, [x(t),y(t)],ti..tf,plotopts,style=point,symbol=circle,symbolsize=5,frames=N,pwrct);Poincare section \342\200\223 brute force version # may be better if more elegant approach complains about maxfun N := 60; # number of periods to plot ni := 30; # starting period of section # add (ni,N) to parameter list in plot title plist := cat(plist,"\n", "ni=",convert(ni,string)," ", "N=",convert(N,string) ) ; # sequence of points consisting of [x,y] pairs at times ni*T to (ni+N)*T; output suppressed sectionpoints := { seq( [ rhs(numsol(ni*T + i*T)[2]) , rhs(numsol(ni*T + i*T)[3]) ], i=0..N-1) }: # same sequence but with x mod 2Pi # sectionpoints := { seq( [ mod2Pi(rhs(numsol(ni*T + i*T)[2])) , rhs(numsol(ni*T + i*T)[3]) ], i=0..N-1) }: pointplot( sectionpoints ,symbol=circle,symbolsize=1, plotopts, pwrct);To study Lyapunov exponents/dependence on initial conditions we need another set of initial conditions and their numerical solution: x1 := x0+0.1; y1 := y0; z1 := z0; ic1 := x(0)=x1,y(0)=y1,z(0)=z1; list of parameter values for title/legend plist1 := cat(plist," ", "x1=",convert(x1,string)," ", "y1=",convert(y1,string) ) ; # options for appearance of plots plotopts1 := 'axes=BOXED,title=plist1,titlefont=[HELVETICA,16]' ; Numerical solution with initial conditions ic1 numsol1 := dsolve( [eq1,eq2,eq3,ic1],[x(t),y(t),z(t)],numeric ); Difference between neighboring solutions # Only way I can so far figure out how to define difference of solutions \342\200\223 some twisted thing with map? # The quotes '' around the definition are to prevent "plot" from attempting to evaluate "rhs" before # plot has told it what t should be. Dnumsol := t -> 'rhs( numsol(t)[2] ) - rhs( numsol1(t)[2] )';Log of difference between neighboring solutions plot( log(abs(Dnumsol(t))) , t=0..15,plotopts1,numpoints=400); Plot both solutions on same graph # clumsy version # numsolX := t -> 'rhs( numsol(t)[2] )'; # numsolX1 := t -> 'rhs( numsol1(t)[2] )'; # plot( {numsolX(t),numsolX1(t)} , t=0..25,plotopts1,numpoints=400); # This is the ticket: pt := odeplot( numsol, [t,x(t)],0..1,plotopts1,color=red,numpoints=400): pt1 := odeplot( numsol1, [t,x(t)],0..1,plotopts1,color=blue,numpoints=400): display({pt,pt1});Plot both solutions on same phase graph # clumsy way: numeric implicit plot # numsolV := t -> 'rhs( numsol(t)[3] )'; # numsolV1 := t -> 'rhs( numsol1(t)[3] )'; # plot( [[numsolX(s),numsolV(s), s=0..25],[numsolX1(t),numsolV1(t), t=0..25]] , plotopts1,numpoints=400); # This is the ticket: ptP := odeplot( numsol, [[x(t),y(t)]],0..1,plotopts1,color=red,numpoints=400): ptP1 := odeplot( numsol1, [[x(t),y(t)]],0..1,plotopts1,color=blue,numpoints=400): display({ptP,ptP1});