(* :Title: SpringSystem *) (* :Author: Christopher Moretti *) (* :Summary: This package contains the command SpringSystem, which creates plots of an arbitrary system of masses and springs for animation. *) (* :History: Version 1.0 by Christopher Moretti, 1997 *) Clear[SpringSystem] SpringSystem::usage = "SpringSystem[parameters,start,finish,frames,maxstep] creates a sequence of \ plots of a 1-dimensional system of springs and masses along with a meter of the \ relative total energy. The plots cover the time from t='start' to t='finish' \ with 'frames' number of plots per second. The list \ 'parameters'={list1,list2...} defines the constants for the system, with \ 'listi' defining the ith mass and spring as follows: {mass,friction \ constant,spring constant,rest position,initial displacement,initial \ velocity}. The option MaxSteps is the maximum number of steps allowed in the \ numerical solution by NDSolve. The option ExternalForce is the force (which \ must be a function of t) which acts on the rightmost mass. Its default value \ is 0." Options[SpringSystem]={MaxSteps->1000,ExternalForce->0} SpringSystem[param_List,start_,finish_,frames_,opts___]:= Block[{m,d,k,l,y0,y1,poslist,vellist,enlist, totmax,totmin,rad,max1,min1,sols,f, enmax,enmeter,y, eqnlis,funclist,a,base,wall,masslist, templist,enframe,n,dt,partrans,j,maxstep,extfor}, (* Parameters: m[i]:ith mass d[i]:coefficient of friction for ith mass k[i]:spring constant for ith spring l[i]:rest position for ith mass y0[i]:initial displacement of ith mass y1[i]:initial velocity of ith mass poslist: matrix of positions for masses; poslist[[i]]=list of positions for ith mass vellist: matrix of velocities for masses (numerical); vellist[[i]]=list of velocities for ith mass enlist: list of total energy (mechanical and potential) for the system. max1: largest position to the right min1: the smallest position rad: radius used for plotting the masses totmax: largest x-value for plotting the masses totmin: the smallest x-value (0 included) for plotting the masses sols: the solution returned by NDSolve for the system of differential equations. f[i][t]: the position of the ith mass at time t enmax: the maximum total energy of the system enmeter: the filled- in rectangle representing the current total energy of the system y[i][t]: displacement of the ith mass from rest position at time t- used in the differential equations eqnlis: list of differential equations funclist: list of the functions y[i] a: iterator base: the "floor" the masses rest on in the plots wall: the wall (at position 0) to which spring 1 is attached masslist: list of the graphics objects representing the masses templist: list of the graphics objects for the springs enframe: the graphics object for the frame of the energy meter n: number of plots dt: time incremement between plots partrans: transpose of the parameter list - used in checking for positive/nonnegative parameters j: number of equations maxstep: maxmimum number of steps for NDSolve extfor: external force on last mass *) (* Parameter Check *) partrans=Transpose[param]; If[MemberQ[Positive[partrans[[1]]],False], Print["Nonpositive mass encountered."];Return[]]; If[MemberQ[Negative[partrans[[2]]],True], Print["Negative friction coefficient detected"];Return[]]; If[MemberQ[Positive[partrans[[3]]],False], Print["Nonpositive spring constant encountered."];Return[]]; (* Define j,dt,n,funclist,m,d,k,l,y0,y1,maxstep,extfor *) j=Length[param]; dt=1/frames //N; n=Max[1,Ceiling[frames(finish-start)]]; funclist=Table[y[i],{i,1,j}]; y[0][t_]:=0;l[0]=0; Do[{m[i],d[i],k[i],l[i],y0[i],y1[i]}= param[[i]],{i,1,j}]; {maxstep,extfor}={MaxSteps,ExternalForce}/.{opts}/ .Options[SpringSystem]; (* Set up equations, solve them, define functions f[i][t] *) eqnlis=Flatten[ Union[{m[j]y[j]''[t]==-d[j]y[j]'[t]- (y[j][t]-y[j-1][t])k[j]+extfor,y[j][0]= =y0[j], y[j]'[0]==y1[j]},Table[ {m[i]y[i]''[t]==-d[i]y[i]'[t]- (y[i][t]-y[i-1][t])k[i]+ (y[i+1][t]-y[i][t])k[i+1], y[i][0]==y0[i],y[i]'[0]==y1[i]},{i,1,j-1}] ]]; sols=NDSolve[eqnlis,funclist,{t,start-1,finish+1}, MaxSteps->maxstep]; Do[f[i][t_]:=Evaluate[y[i][t]/. sols],{i,1,j}]; (* Create matrixes/lists for positions, velocities, energy *) poslist=Table[Flatten[Table[f[i][start+a dt]+l[i],{a,0,n}]], {i,1,j}]; vellist=Table[Flatten[Table[ (f[i][start+a dt+.001]-f[i][start+a dt-.001])/.002,{a,0,n} ]],{i,1,j}]; enlist=k[1] (poslist[[1]]-l[1])^2/2+m[1]vellist[[1]]^2/2+ Sum[k[i](poslist[[i]]-l[i]-poslist[[i-1]]+l[i-1])^2+ m[i] vellist[[i]]^2,{i,2,j}]/2; (* Define enmax, max1, min1, totmax, totmin, rad *) enmax=Max[Flatten[enlist]]; max1=Max[Abs[Flatten[poslist]]]; min1=Min[Flatten[poslist],{0}]; rad=max1/20;totmax=max1+rad;totmin=min1-2rad; (* Create the graphics objects base,wall, enframe - constant from plot to plot *) base=Polygon[{{totmin,-rad},{totmax,-rad},{totmax,- 7rad/6}, {totmin,-7/6 rad}}]; wall=Polygon[{{totmin,-7 rad/6},{totmin,rad}, {totmin-rad/6,rad},{totmin-rad/6,-7rad/6}}]; enframe=Graphics[Line[{{totmax+22 rad/8,-rad}, {totmax+27 rad/8,-rad},{totmax+27 rad/8,rad}, {totmax+22rad/8,rad},{totmax+22 rad/8,- rad}}]]; (* Do loop to create plots - creates masslist,templist, and enmeter and then combines them all in a single plot with Show *) Do[ masslist= Table[ Graphics[Disk[{poslist[[i,a]],0},rad]],{i,1,j}]; templist = Union[ {Plot[ rad Sin[(t-totmin)10Pi/(poslist[[1,a]] -totmin)],{t,totmin,poslist[[1,a]]}, DisplayFunction->Identity]}, Table[Plot[ rad Sin[(t+poslist[[i,a]])10Pi/(poslist[[i+1,a]] -poslist[[i,a]])],{t,poslist[[i,a]], poslist[[i+ 1,a]]},DisplayFunction->Identity], {i,1,j-1}]]; enmeter=Graphics[Rectangle[{totmax+22rad/8,-rad}, {totmax +27 rad/8,2 rad enlist[[a]]/enmax- rad}]]; Show[masslist,Graphics[base],Graphics[wall],templist,enmeter,enframe, AspectRatio->Automatic, PlotRange->{{totmin-rad,totmax+4rad},{-1.1rad,rad}} ]; ,{a,1,n}]; ]; (* Examples SpringSystem[{{1,1,4,0,2,6}},0,10,1, ShowMeter->True,ShowRestPositions->True] SpringSystem[{{1,0,2,3,1,0},{1,0,4,7,-1,0}},0,10,5] SpringSystem[{{1,0,2,3,1,0},{1,0,4,7,-1,0},{1,0,1,10,2,0}}, 0,20,10,ExternalForce->Sin[t],ShowRestPositions->True] SpringSystem[{{1,0,4,0,1,0}},0,50,5, ExternalForce->Sin[2t]/4,ShowRestPositions->True] *)