Spring simulation

Oscillating springs, damped springs, and driven springs

- Marcus McGuff, Austin Community College

If you want to enter your own driving, damping, or spring force function and you are running a full version of Mathematica, you can click on the little triangle next to each of those to open up a box you can type directly into (either using standard mathematical form or Mathematica notation).  You can set initial conditions using either the y(0)/v(0) sliders or by clicking your mouse on the y/y' phase plane.  If you want to "zoom in" to look more carefully at the transient response, click on "shorter time"; if you want to "zoom out" to examine more of the long-term behavior, click on "longer time".  The driving force must be only a function of time, t; any other "variables" that are in that formula will be treated as parameters.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

plotRangeLimits[plotGraph_,component_:2]:=If[Length[Cases[plotGraph,GraphicsComplex[pt_,stuff__]->pt,Infinity]]>0,Through[{(Min[#]&),(Max[#]&)}[Flatten[Cases[plotGraph,GraphicsComplex[pt_,stuff__]->pt,Infinity],1][[All,component]]]],Through[{(Min[#]&),(Max[#]&)}[Flatten[Cases[plotGraph,Line[pt_]->pt,Infinity],1][[All,component]]]]]

fixDomainRange[dRange_]:=Module[{domainRange=dRange},
Table[{
If[domainRange[[component,2]]<0,domainRange[[component,2]]=0];
If[domainRange[[component,1]]>0,domainRange[[component,1]]=0];
If[domainRange[[component,1]]==domainRange[[component,2]],domainRange[[component]]={-1,1}];},{component,1,2}];
1.1*domainRange]
(* This just makes sure the domain and range always include 0 and fixes the possibility where min=max value*)

ticks[min_,max_]:=Table[i,{i,Floor[max],Ceiling[min],Round[(Ceiling[min]-Floor[max])/4]}]

springFunctions={
{ToBoxes[k y]},
{ToBoxes[( k(1-(E^(-ϵ t)))+k0 (E^(-ϵ t)))y]},
{ToBoxes[ k y-ϵ (y^3)]}
};

dampingFunctions={
{ToBoxes[b y']}
};

drivingFunctions:={
{ToBoxes[A Sin[ω t]]},
{ToBoxes[A Sin[ω t-ϕ]]},
{ToBoxes[A (E^(-k t))]}(*,
{ToBoxes[A DiracDelta[t]]}*)
};

lengthSpring=6;ωSpring=16;maxStretch=0.9;

spring={First[ParametricPlot3D[{Sin[ωSpring u],Cos[ωSpring u],u },{u,0,lengthSpring},PlotStyle->{Orange,Opacity[1],Thickness[0.015]}]]};

Manipulate[
If[var===Null,Clear[var]];

(* Initialize list of functions *)
If[springText===Null,springText=springFunctions[[1,1]]];
If[drivingText===Null,drivingText=drivingFunctions[[1,1]]];
If[dampingText===Null,dampingText=dampingFunctions[[1,1]]];
If[scaledStretch===Null,scaledStretch=maxStretch];

(*springODEq=If[dampingYN,y''[t]+(springEq/.y'->yTMP'[t]/.y->yTMP[t]/.yTMP->y/.springParams)+(dampingEq/.y'->yTMP'[t]/.y->yTMP[t]/.yTMP->y/.dampingParams),y''[t]+(springEq/.y'->yTMP'[t]/.y->yTMP[t]/.yTMP->y/.springParams)]-If[drivenYN,(drivingEq/.y'->yTMP'[t]/.y->yTMP[t]/.yTMP->y/.drivingParams),0];*)

(*Refresh[(*Print["Refresh loop"]*)*)
springEq=ToExpression[springText,TraditionalForm];
springParams=SortBy[Thread[Rule[#,Replace[#,Join[springParams,{_->1.}],1]]],First]&[Cases[Level[{springEq},-1],s_Symbol/;And[FreeQ[{t,y,y'},s],FreeQ[Attributes[s],Locked|Constant]]]];

drivingEq=ToExpression[drivingText,TraditionalForm];
drivingParams=SortBy[Thread[Rule[#,Replace[#,Join[drivingParams,{_->1.}],1]]],First]&[Cases[Level[{drivingEq},-1],s_Symbol/;And[FreeQ[{t},s],FreeQ[Attributes[s],Locked|Constant]]]];
drivingPos=Integrate[Integrate[drivingEq/.t->tmp,tmp],tmp]/.tmp->t;

dampingEq=ToExpression[dampingText,TraditionalForm];
dampingParams=SortBy[Thread[Rule[#,Replace[#,Join[dampingParams,{_->1.}],1]]],First]&[Cases[Level[{dampingEq},-1],s_Symbol/;And[FreeQ[{t,y,y'},s],FreeQ[Attributes[s],Locked|Constant]]]];

ySol[tFake_]=First[yTmp[tFake]/.
NDSolve[{
((If[dampingYN,
 m*yTmp''[tFake]+(springEq/.y'->yTmp'[tFake]/.y->yTmp[tFake]/.springParams)+(dampingEq/.y'->yTmp'[tFake]/.y->yTmp[tFake]/.dampingParams),
m*yTmp''[tFake]+(springEq/.y'->yTmp'[tFake]/.y->yTmp[tFake]/.springParams)]-If[drivenYN,(drivingEq/.y'->yTmp'[tFake]/.y->yTmp[tFake]/.drivingParams),0])/.t->tFake)
==0,
yTmp[0]==y0,yTmp'[0]==v0},
yTmp[tFake],{tFake,0,tLimit}]];

ySolD[tFake_]=ySol'[tFake];

var["yDR"]={Flatten@InterpolatingFunctionDomain[Head@ySol[tmp]],{Min[InterpolatingFunctionValuesOnGrid[Head@ySol[tmp]]],Max[InterpolatingFunctionValuesOnGrid[Head@ySol[tmp]]]}};
yDomainRange=fixDomainRange[var["yDR"]];

(*var["dDR"]={{0,tLimit},{NMinValue[{-ySol[t],0≤t≤tLimit},t],NMaxValue[{-ySol[t],0≤t≤tLimit},t]}};
derivDomainRange=fixDomainRange[var["dDR"]];

var["pDR"]={yDomainRange[[2]],derivDomainRange[[2]]};
phaseDomainRange=fixDomainRange[var["pDR"]];
*)
var["derivPlot"]=Plot[ySolD[tmp],{tmp,0,tLimit},PlotRange->All];
var["dDR"]={{0,tLimit},plotRangeLimits[var["derivPlot"]]};
derivDomainRange=fixDomainRange[var["dDR"]];

var["phasePlot"]=ParametricPlot[{ySol[tmp],ySolD[tmp]},{tmp,0,tLimit}];
var["pDR"]={plotRangeLimits[var["phasePlot"],1],plotRangeLimits[var["phasePlot"],2]};
phaseDomainRange=fixDomainRange[var["pDR"]];

If[drivenYN,
var["forcingPlot"]=Plot[(drivingEq/.drivingParams/.t->tmp),{tmp,0,tLimit},PlotRange->All];
var["fDR"]={{0,tLimit},plotRangeLimits[var["forcingPlot"]]};
forcingDomainRange=fixDomainRange[var["fDR"]];
var["forcingPosPlot"]=Plot[(drivingPos/.drivingParams/.t->tmp),{tmp,0,tLimit},PlotRange->All];
var["fpDR"]={{0,tLimit},plotRangeLimits[var["forcingPosPlot"]]};
forcingPosDomainRange=fixDomainRange[var["fpDR"]]
];
(*,
TrackedSymbols:>{dampingYN,drivenYN,springText, drivingText,dampingText,y0,v0,tLimit}];*)

With[{
yMax=Max[Abs[yDomainRange[[2]]]],
fMax=Max[Abs[forcingPosDomainRange[[2]]]]},
scaledStretch=(maxStretch*If[drivenYN,
yMax/(yMax+fMax),1])];

Grid[{{Text@Style[Row[{"m y''",
If[dampingYN,
Row[{" + ",ToExpression[dampingText,StandardForm,TraditionalForm]," + "}],
" + "],
ToExpression[springText,StandardForm,TraditionalForm]," = ",If[drivenYN,ToExpression[drivingText,StandardForm,TraditionalForm],0]}]],SpanFromLeft},

With[{
(* These control the position of the end of the spring and the piston when forcing *)
endSpring=scaledStretch*((-ySol[time]/Max[Abs[yDomainRange[[2]]]])+1),
endForcing=lengthSpring*(maxStretch-scaledStretch)*((-drivingPos/.drivingParams/.t->time)/(Max[Abs[forcingPosDomainRange[[2]]]])+1)},

{Dynamic@Column[{

(* Plot forcing, y[t], y'[t] on left *)
If[drivenYN,
Plot[{ySol[t],(drivingEq/.drivingParams)},{t,0,Max[time,0.0001]},AxesLabel->{"t",Column[{Style["y(t)",Red],Style["F(t)",Blue]}]},PlotRange->{forcingDomainRange[[1]],{Min[forcingDomainRange[[2]],yDomainRange[[2]]],Max[forcingDomainRange[[2]],yDomainRange[[2]]]}},PlotStyle->{Directive[Red,Thickness[0.01]],Directive[Blue,Thickness[0.01],Dashing[Small]]},Epilog->{Red,PointSize[Large],Point[{time,ySol[time]}],Blue,Point[{time,(drivingEq/.drivingParams/.t->time)}]}],Plot[ySol[t],{t,0,Max[time,0.0001]},AxesLabel->{"t","y(t)"},PlotRange->yDomainRange,PlotStyle->Directive[Red,Thickness[0.01]],Epilog->{Red,PointSize[Large],Point[{time,ySol[time]}]}]],,
Plot[ySolD[t],{t,0,Max[time,0.0001]},AxesLabel->{"t","y'(t)"},PlotRange->derivDomainRange,PlotStyle->Directive[Green,Thickness[0.01]],Epilog->{Green,PointSize[Large],Point[{time,ySolD[time]}]}], ,
If[drivenYN,
Plot[ySol[t]+(drivingPos/.drivingParams),{t,0,Max[time,0.0001]},AxesLabel->{"t",Column[{"Position:","Spring+piston"},Alignment->Center]},PlotRange->{forcingDomainRange[[1]],forcingDomainRange[[2]]+yDomainRange[[2]]},PlotStyle->Directive[Gray,Thickness[0.01]],Epilog->{Gray,PointSize[Large],Point[{time,ySol[time]+(drivingPos/.drivingParams/.t->time)}]}],""]
}],,

Dynamic@Column[{Graphics3D[{
(* Despite picture labels, +z (y) is downwards *)
(* Draw piston if driven, rectangular plate if not *)
If[drivenYN,
{Cylinder[{{0,0,0},{0,0,endForcing}},0.3],Cylinder[{{0,0,endForcing},{0,0,endForcing+0.1}},1.3]
},
Polygon[2{{-1,-1,0},{-1,1,0},{1,1,0},{1,-1,0}}]],

(* If driven, make spring start at end of piston, then draw spring *)
If[drivenYN,
Translate[Scale[spring,{1,1,endSpring},{0,0,0}],{0,0,endForcing}],
Scale[spring,{1,1,endSpring},{0,0,0}]],

(* Draw/label eq point for spring relative to piston *)
Arrowheads[Small],
Blue,
If[drivenYN,Arrow[{{-3.5,0,If[drivenYN,2+endForcing,0]},{-3.5,0,If[drivenYN,2+endForcing-2(drivingEq/.drivingParams/.t->time)/Max[Abs[forcingDomainRange[[2]]]],0]}}],{}],
If[drivenYN,Text["(Overscript[F, ->])",{-4,0,If[drivenYN,1.5+endForcing,0]}],{}],
Arrowheads[Medium],
Gray,
Arrow[{{-1.2,0,If[drivenYN,endForcing,0]},{-3,0,If[drivenYN,endForcing,0]},{-3,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]},{-1.2,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]}}],
Text["Eq.",{-4.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]}],
Arrow[{{-5.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]+1.5},
{-5.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]+2.5}}],
Text["-y",{-5.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]+1}],
Arrow[{{-5.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]-1.5},
{-5.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]-2.5}}],
Text["+y",{-5.5,0,lengthSpring*scaledStretch+If[drivenYN,endForcing,0]-1}],

(* Draw mass *)
LightGray,
Cylinder[{{0,0,endSpring*lengthSpring+If[drivenYN,endForcing,0]},{0,0,endSpring*lengthSpring+If[drivenYN,endForcing,0]+(m/10)1.5}},1.2],
Black,
Text["m",{0,-1.2,endSpring*lengthSpring+If[drivenYN,endForcing,0]+0.35}],
Arrowheads[Small],
Green,
Arrow[{{2,0,lengthSpring*endSpring+If[drivenYN,endForcing,0]},
{2,0,lengthSpring*endSpring+If[drivenYN,endForcing,0]-2ySolD[time]/Max[Abs[derivDomainRange[[2]]]]}}],
Text["(Overscript[v, ->])",{2.5,0,lengthSpring*endSpring+If[drivenYN,endForcing,0]}],


(* If damped, draw tube/fluid *)
If[dampingYN,
{Opacity[0.2],
Cylinder[{{0,0,0.5},{0,0,1}},1.5],Opacity[0.3],
Cylinder[{{0,0,1},{0,0,(maxStretch+1)*lengthSpring+0.5}},1.5],Opacity[1]},
{}],

(* Draw "ground" *)
Polygon[{{-2,-2,(maxStretch+1)*lengthSpring+0.5},{-2,2,(maxStretch+1)*lengthSpring+0.5},{2,2,(maxStretch+1)*lengthSpring+0.5},{2,-2,(maxStretch+1)*lengthSpring+0.5}}]

},ImageSize->120,Boxed->False,ViewPoint->Front,ViewVertical->{0,0,-1}],


EventHandler[
ParametricPlot[{ySol[t],ySolD[t]},{t,0,Max[0.001,time]},PlotRange->phaseDomainRange,AxesLabel->{"y","y'"},ImageSize->200,PlotStyle->Thickness[0.008],TicksStyle->Opacity[0.8],Epilog->{Yellow,PointSize[Large],Point[{y0,v0}],Red,Point[{ySol[time],ySol'[time]}]}],
{"MouseClicked":>(
{y0,v0}=MousePosition["Graphics"];
time=0;)
}]
}]
}]
}],

{{tLimit, 40}, ControlType -> None}, 

Row[{Column[{
Row[{"         Time: ",Dynamic@time}],Control@{{time, 0,"Time:"},0,Dynamic@tLimit, 0.1,Appearance->{"Labeled","Open"},AnimationRepetitions->1,DefaultDuration->20,AnimationDirection->Forward,AnimationRunning->False,ControlType->Animator,AppearanceElements->{"ProgressSlider","ResetButton","StepLeftButton","PlayPauseButton","StepRightButton","FasterSlowerButtons"}},
Row[{"        ",
Button["<< Shorter time",tLimit=tLimit/2;time=0,BaseStyle->Small,ImageSize->All],
"                ",
Button["Longer time >>",tLimit=tLimit*2,BaseStyle->Small,ImageSize->All]}]}],
"   ",
Column[{Control@{{m, 5,"mass ="},0.5,10, Appearance->{"Labeled","Closed"},ImageSize->Small},Control@{{y0, 1,"y(0) = "},-5,5, Appearance->{"Labeled","Closed"},ImageSize->Small},
Control@{{v0, 0,"y'(0) = "},-5,5, Appearance->{"Labeled","Closed"},ImageSize->Small}}]
}],

{{springParams,{}},ControlType->None},
{{dampingParams,{}},ControlType->None},
{{drivingParams,{}},ControlType->None},

{springEq,ControlType->None},
{dampingEq,ControlType->None},
{drivingEq,ControlType->None},
{drivingPos,ControlType->None},

{{yDomainRange,{{0,tLimit},{-1,1}}},ControlType->None},
{{forcingDomainRange,{-1,1}},ControlType->None},
{{derivDomainRange,{-1,1}},ControlType->None},
{{phaseDomainRange,{-1,1}},ControlType->None},
{{forcingPosDomainRange,{-1,1}},ControlType->None},

Item[
Column[{
Column[{
OpenerView[{Control[{{springText,,"Spring force ="},(#->ToExpression[#,StandardForm,TraditionalForm])&/@springFunctions[[All,1]],ControlType->PopupMenu}],
Dynamic@If[$LicenseType≠"Player",
Control[{{springText,,"Spring force ="},TraditionalForm[InputField[Dynamic[springText],Boxes]]&}],
Text@Style[
"To enter your own equations, you must use a full version of Mathematica",Italic]]}],
Column[{
"Parameters",Dynamic[Grid[MapIndexed[{springParams[[First@#2,1]],Manipulator[Dynamic[springParams[[First@#2,2]]],{0,10,0.01},Appearance->{"Labeled","Closed"},ImageSize->Medium]}&,springParams]]],}]
}],

Dynamic@If[dampingYN,Column[{Control@{{dampingYN,False,"Damping?"},{False,True}},
OpenerView[{Control[{{dampingText,,"Damping ="},(#->ToExpression[#,StandardForm,TraditionalForm])&/@dampingFunctions[[All,1]],ControlType->PopupMenu}],
Dynamic@If[$LicenseType≠"Player",
Control[{{dampingText,,"Damping ="},TraditionalForm[InputField[Dynamic[dampingText],Boxes]]&}],
Text@Style[
"To enter your own equations, you must use a full version of Mathematica",Italic]]}],
Column[{
"Parameters",Dynamic[Grid[MapIndexed[{dampingParams[[First@#2,1]],Manipulator[Dynamic[dampingParams[[First@#2,2]]],{0,10,0.01},Appearance->{"Labeled","Closed"},ImageSize->Medium]}&,dampingParams]]],}]
}],
Control@{{dampingYN,False,"Damping?"},{False,True}}],

Dynamic@If[drivenYN,Column[{Control@{{drivenYN,False,"Driving force?"},{False,True}},
OpenerView[{Control[{{drivingText,,"F(t) ="},(#->ToExpression[#,StandardForm,TraditionalForm])&/@drivingFunctions[[All,1]],ControlType->PopupMenu}],
Dynamic@If[$LicenseType≠"Player",
Control[{{drivingText,,"F(t) ="},TraditionalForm[InputField[Dynamic[drivingText],Boxes]]&}],
Text@Style[
"To enter your own equations, you must use a full version of Mathematica",Italic]]}],
Column[{
"Parameters",Dynamic[Grid[MapIndexed[{drivingParams[[First@#2,1]],Manipulator[Dynamic[drivingParams[[First@#2,2]]],{0,10,0.01},Appearance->{"Labeled","Closed"},ImageSize->Medium]}&,drivingParams]]],}]
}],
Control@{{drivenYN,False,"Driving force?"},{False,True}}]},Dividers->{False,{True,True,True,True}},Alignment->Left],ControlPlacement->Right],

Delimiter,

{endSpring,ControlType->None},
{endForcing,ControlType->None},
{scaledStretch,ControlType->None},

(*temporary holding variable*)
{var,ControlType->None},

TrackedSymbols:>{dampingYN,drivenYN,springText, drivingText,dampingText,springParams,dampingParams,drivingParams,y0,v0,tLimit,time,m},

(*Initialization:>(
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
plotRangeLimits[plotGraph_,component_:2]:=If[Length[Cases[plotGraph,GraphicsComplex[pt_,stuff__]->pt,Infinity]]>0,Through[{(Min[#]&),(Max[#]&)}[Flatten[Cases[plotGraph,GraphicsComplex[pt_,stuff__]->pt,Infinity],1][[All,component]]]],Through[{(Min[#]&),(Max[#]&)}[Flatten[Cases[plotGraph,Line[pt_]->pt,Infinity],1][[All,component]]]]];
fixDomainRange[dRange_]:=Module[{domainRange=dRange},
Table[{
If[domainRange[[component,2]]<0,domainRange[[component,2]]=0];
If[domainRange[[component,1]]>0,domainRange[[component,1]]=0];
If[domainRange[[component,1]]==domainRange[[component,2]],domainRange[[component]]={-1,1}];},{component,1,2}];
1.1*domainRange]
(* This just makes sure the domain and range always include 0 and fixes the possibility where min=max value*);
springFunctions={
{ToBoxes[k y]},
{ToBoxes[( k(1-(E^(-ϵ t)))+k0 (E^(-ϵ t)))y]},
{ToBoxes[ k y-ϵ (y^3)]}
};
dampingFunctions={
{ToBoxes[b y']}
};
drivingFunctions:={
{ToBoxes[A Sin[ω t]]},
{ToBoxes[A Sin[ω t-ϕ]]},
{ToBoxes[A (E^(-k t))]}(*,
{ToBoxes[A DiracDelta[t]]}*)
};
lengthSpring=6;ωSpring=18;maxStretch=0.8;
spring={First[ParametricPlot3D[{Sin[ωSpring u],Cos[ωSpring u],u },{u,0,lengthSpring},PlotStyle->{Orange,Opacity[1],Thick}]]}
)*)

SaveDefinitions->True]