Nest numerical algorithms

Harald Biller, Continental Automotive Technologies

With[
{g=9.81,m=0.75,r=0.02,l=0.05,J=(10^(-4)),δ=0.001},
With[
{
sol=
Block[
{ϕ,T},
{
NDSolveValue[
{
ϕ''[t]==-((m (r^2) Sin[ϕ[t]] (g+((((r^2) Cos[ϕ[t]])/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))])))+(((r^4) ((Sin[ϕ[t]])^2))/((l^3) (((1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))))^(3/2))))) ((ϕ'[t])^2)))/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))]) (J+((m (r^4) ((Sin[ϕ[t]])^2))/((l^2) (1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2)))))))),
ϕ[0]==(π/2),ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,T=t;"StopIntegration"]
},
ϕ,{t,0,100}
],
T
}
]
},
Manipulate[
Module[
{ϕ,T,h,fixed,moving},
{ϕ,T}=sol;
ϕ=ϕ[τ T];
h=l(1-(Sqrt[1-2((((r/l)))^2)(1-Cos[ϕ])]));
{fixed,moving}=Transpose @ Table[{{0,-r,l},{r Sin[ϕ],-r Cos[ϕ],h}}.RotationMatrix[ψ,{0,0,1}],{ψ,0,((4π)/3),((2π)/3)}];
Graphics3D[
{PointSize[0.02],Point[{0,0,0}],Thickness[0.005],EdgeForm[Thickness[0.005]],Line[Table[{r Cos[ψ],r Sin[ψ],0},{ψ,0,2π,(π/20)}]],
If[show,{{Thickness[0.0025],Line[{{0,-r,0},{0,-r,l}}]},Line[{(*{fixed[[1]],{0,-r,h},moving[[1]]},*){{0,0,h},{0,-r,h}},{{0,0,h},(({0,-r,h}+moving[[1]])/2)}(*,Table[{r Sin[ψ],-r Cos[ψ],h},{ψ,0,ϕ,(ϕ/20)}]*)}],Opacity[0.5],Polygon[{fixed[[1]],{0,-r,h},moving[[1]]}],Polygon[Append[Table[{r Sin[ψ],-r Cos[ψ],h},{ψ,0,ϕ,(ϕ/20)}],{0,0,h}]]},{}],Red,Tube[Map[{#,{0,0,h}}&,moving],δ],Blue,MapThread[Tube[{#1,#2},δ]&,{fixed,moving}]},
PlotRange->{{-r-2δ,r+2δ},{-r-2δ,r+2δ},{-2δ,l+2δ}},ViewPoint->{1.75,-2.4,1.6},
ImageSize->400
]
],
(*{{τ,0.25},0,1,0.025,Appearance->"Labeled"},*)
{τ,0,1,ControlType->Animator},
{{show,False,"show auxiliary lines"},ControlType->None},
LabelStyle->24
]
]
]

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Rotation pendulum

Observe a rotation pendulum to obtain a body’s moment of inertia.

With[
{g=9.81,m=0.75,r=0.02,l=0.05,J=(10^(-4)),δ=0.001},
With[
{
sol=
Block[
{ϕ,T},
{
NDSolveValue[
{
ϕ''[t]==-((m (r^2) Sin[ϕ[t]] (g+((((r^2) Cos[ϕ[t]])/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))])))+(((r^4) ((Sin[ϕ[t]])^2))/((l^3) (((1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))))^(3/2))))) ((ϕ'[t])^2)))/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))]) (J+((m (r^4) ((Sin[ϕ[t]])^2))/((l^2) (1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2)))))))),
ϕ[0]==(π/2),ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,T=t;"StopIntegration"]
},
ϕ,{t,0,100}
],
T
}
]
},
Manipulate[
Module[
{ϕ,T,h,fixed,moving},
{ϕ,T}=sol;
ϕ=ϕ[τ T];
h=l(1-(Sqrt[1-2((((r/l)))^2)(1-Cos[ϕ])]));
{fixed,moving}=Transpose @ Table[{{0,-r,l},{r Sin[ϕ],-r Cos[ϕ],h}}.RotationMatrix[ψ,{0,0,1}],{ψ,0,((4π)/3),((2π)/3)}];
Graphics3D[
{PointSize[0.02],Point[{0,0,0}],Thickness[0.005],EdgeForm[Thickness[0.005]],Line[Table[{r Cos[ψ],r Sin[ψ],0},{ψ,0,2π,(π/20)}]],
If[show,{{Thickness[0.0025],Line[{{0,-r,0},{0,-r,l}}]},Line[{(*{fixed[[1]],{0,-r,h},moving[[1]]},*){{0,0,h},{0,-r,h}},{{0,0,h},(({0,-r,h}+moving[[1]])/2)}(*,Table[{r Sin[ψ],-r Cos[ψ],h},{ψ,0,ϕ,(ϕ/20)}]*)}],Opacity[0.5],Polygon[{fixed[[1]],{0,-r,h},moving[[1]]}],Polygon[Append[Table[{r Sin[ψ],-r Cos[ψ],h},{ψ,0,ϕ,(ϕ/20)}],{0,0,h}]]},{}],Red,Tube[Map[{#,{0,0,h}}&,moving],δ],Blue,MapThread[Tube[{#1,#2},δ]&,{fixed,moving}]},
PlotRange->{{-r-2δ,r+2δ},{-r-2δ,r+2δ},{-2δ,l+2δ}},ViewPoint->{1.75,-2.4,1.6}
]
],
{{τ,0.25},0,1,0.025,Appearance->"Labeled"},
{{show,False,"show auxiliary lines"},{True,False}},
LabelStyle->24
]
]
]

Clear[m,J,l,r,g,h,ϕ,f,t,T,ϕInit]
constrainedHeight=ϕ|->l(1-(Sqrt[1-2((((r/l)))^2)(1-Cos[ϕ])]));

Equation of motion

defaultParameters={m->0.5,l->1,r->0.02,g->9.81};

energy=(1/2)J ((ϕ'[t])^2)+(1/2)m ((h'[t])^2)+m g h[t]/.h->(t|->f[ϕ[t]]);

conservation of energy:

equation=Simplify @ Solve[0==D[energy,t],ϕ''[t]]

rhs=ϕ''[t]/.equation[[1]]/.f->constrainedHeight

For l->∞, we obtain the mathematical pendulum, ϕ''=-((Subscript[ω, 0])^2)sin(ϕ):

rhsLargeL=Normal @ Echo @ Series[rhs,{l,∞,1}];

Period

Solve equation of motion and detect period:

Block[
{ϕ,t,T,tEnd=1.9,style=24,size=800,ar=(1/4),pad={{50,150},{10,10}}},
ϕ=NDSolveValue[
{ϕ''[t]==-((m (r^2) Sin[ϕ[t]] (g+((((r^2) Cos[ϕ[t]])/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))])))+(((r^4) ((Sin[ϕ[t]])^2))/((l^3) (((1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))))^(3/2))))) ((ϕ'[t])^2)))/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))]) (J+((m (r^4) ((Sin[ϕ[t]])^2))/((l^2) (1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))))))))/.{m->0.5`,l->1,r->0.02`,g->9.81`}/.J->(10^(-4)),ϕ[0]==90 °,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,T=t]},
ϕ,{t,0,tEnd}
];
Column[
{
Plot[ϕ[t]/°,{t,0,tEnd},PlotStyle->Thick,(*PlotTheme->"Detailed",PlotLegends->None,*)BaseStyle->style,ImageSize->size,ImagePadding->pad,AspectRatio->ar,PlotLabel->"angle ϕ(t) [°]",PlotRangePadding->{None,Automatic},Prolog->{Thick,Red,Line[{{T,-180},{T,180}}]},AxesLabel->{"time [s]"}],
Plot[ϕ'[t]/°,{t,0,tEnd},PlotStyle->Thick,(*PlotTheme->"Detailed",PlotLegends->None,*)BaseStyle->style,ImageSize->size,ImagePadding->pad,AspectRatio->ar,PlotLabel->"speed ϕ'(t) [°/s]",PlotRangePadding->{None,Automatic},Prolog->{Thick,Red,Line[{{T,-1000},{T,1000}}]},AxesLabel->{"time [s]"}]
}
]
]

NDSolve[
{ϕ''[t]==rhs/.defaultParameters/.J->(10^(-4)),ϕ[0]==90 °,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,Print[t];"StopIntegration"]},
ϕ,{t,0,20}
]

Clear[period]
period[J1_?NumericQ,ϕInit_?NumericQ,param_]:=
Block[
{T},
NDSolve[
{ϕ''[t]==rhs/.param/.J->J1,ϕ[0]==ϕInit,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,T=t;"StopIntegration"]},
ϕ,{t,0,100}
];
T
]

Note that evaluation is deferred until function arguments are numeric.

period[(10^(-4)),90°,defaultParameters]

Plot[period[(10^(-4)),ψ °,defaultParameters],{ψ,0.1,175},]

Get moment of inertia from period

Plot[period[(10^(-4))J,90°,defaultParameters],{J,1,2},]

Block[
{T1=2,J1},
FindRoot[
period[(10^(-4))J1,90°,defaultParameters]==T1,
{J1,1,2},
AccuracyGoal->6,PrecisionGoal->6
]
]

NDSolve is nested within FindRoot.

A wrapper function avoids premature evaluation.

Relax outer tolerances because of inner numerical noise.

Other examples: ConvexOptimization nested within FindMinimum, ...

https://www.wolfram.com/wolfram-u/courses/mathematics/advanced-numerical-computing-math026/

Deploy through Manipulate

With[
{style={FontSize->48,FontFamily->"Calibri"},tMax=100,g=9.81},
Manipulate[
Module[
{ϕ,t,period,J,m=(10^(-3))m0,r=(10^(-3))r0,l=(10^(-3))l0},
period[J_?NumericQ]:=
Module[
{T},
NDSolve[
{ϕ''[t]==-((m (r^2) Sin[ϕ[t]] (g+((((r^2) Cos[ϕ[t]])/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))])))+(((r^4) ((Sin[ϕ[t]])^2))/((l^3) (((1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))))^(3/2))))) ((ϕ'[t])^2)))/(l (Sqrt[1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2))]) (J+((m (r^4) ((Sin[ϕ[t]])^2))/((l^2) (1-((2 (r^2) (1-Cos[ϕ[t]]))/(l^2)))))))),ϕ[0]==ϕInit °,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,T=t;"StopIntegration"]},
ϕ,{t,0,100}
];
T
];
Row[{"moment of inertia: ",Round[100J,0.1]," kg (mm^2)"},BaseStyle->style]/.FindRoot[period[(10^(-4))J]==T,{J,0.1,2,(10^(-6)),(10^3)},AccuracyGoal->6,PrecisionGoal->6]
],
{{ϕInit,90,"amplitude [°]"},5,175,5,Appearance->"Labeled"},
{{T,2.,"period [s]"},0.1,5.,0.1,Appearance->"Labeled"},
Delimiter,
{{m0,500,"mass [g]"},25,1000,25,Appearance->"Labeled"},
{{r0,20,"radius [mm]"},1,50,1,Appearance->"Labeled"},
{{l0,1000,"wire length [mm]"},100,2500,50,Appearance->"Labeled"},
LabelStyle->style
]
]

The end.

Appendix 1: long wire approximation

Clear[periodLargeL]
periodLargeL[J1_?NumericQ,ϕInit_?NumericQ,param_]:=
Block[
{T},
NDSolve[
{ϕ''[t]==rhsLargeL/.param/.J->J1,ϕ[0]==ϕInit,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,T=t;"StopIntegration"]},
ϕ,{t,0,100}
];
T
]

Block[{T1=2,J1},FindRoot[periodLargeL[(10^(-4))J1,90°,defaultParameters]==T1,{J1,1,2},AccuracyGoal->6,PrecisionGoal->6]]

Actually, a mathematical pendulum can be solved in terms of an elliptic function:

(10^4)((g m (r^2))/l)((((T/(4EllipticK[((Sin[((ϕInit)/2)])^2)]))))^2)/.defaultParameters/.{T->2,ϕInit->90°}

Previous result without the long wire approximation:

Block[{T1=2,J1},FindRoot[period[(10^(-4))J1,90°,defaultParameters]==T1,{J1,1,2},AccuracyGoal->6,PrecisionGoal->6]]

But replacing Sin[ϕ]->ϕ, we get:

(10^4)((g m (r^2))/l)((((T/(2π))))^2)/.defaultParameters/.T->2

Appendix 2: calculate the period as an integral

Block[
{J,ϕInit,tϕ,tϕLargeL,TLargeL},
tϕ=Echo[(1/(ϕ'[t]))/.Last @ Solve[(energy/.ϕ[t]->ϕ)==(energy/.{ϕ[t]->ϕInit,ϕ'[t]->0})/.f->constrainedHeight,ϕ'[t]],"integrand:"];
Assuming[
g>0∧l>0∧m>0∧r>0∧J>0∧0<ϕInit<π,
tϕLargeL=Simplify @ Normal @ Series[tϕ,{l,∞,1}];
TLargeL=Echo[Simplify @ Integrate[4tϕLargeL,{ϕ,0,ϕInit}],"long wire approximation of period:"]
];
J=(10^(-4));ϕInit=90°;tϕ=tϕ/.defaultParameters;
Echo[TLargeL/.defaultParameters,"numerical value of period with long wire approximation:"];
Echo[4NIntegrate[tϕ,{ϕ,0,ϕInit}],"numerical value of period without long wire approximation:"];
]

Compare this to the values based on NDSolve and WhenEvent:

{periodLargeL[(10^(-4)),90°,defaultParameters],period[(10^(-4)),90°,defaultParameters]}

Appendix 3: ParametricNDSolve

WhenEvent works with ParametricNDSolve. Note that this function caches results. We Sow the period, and this happens only when a parameter combination is used for the first time, unless we switch off caching by Method -> {“ParametricCaching” -> None}.

Alternatively, we could replace Sow by a suitable mechanism from the package DifferentialEquations`InterpolatingFunctionAnatomy`, or get the final time from the InterpolatingFunction ϕ by ϕ[“Grid”][[-1]].

My first experiment used every parameter at most once, so the problem did not occur:

Block[
{J,ϕ1,period,ϕInit=90°,param=defaultParameters,T=2},
ϕ1=
ParametricNDSolveValue[
{ϕ''[t]==rhs/.param,ϕ[0]==ϕInit,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,Sow[t];"StopIntegration"]},
ϕ,{t,0,20},{J}
];
period[J_?NumericQ]:=Reap[ϕ1[J]][[2,1,1]];
FindRoot[T==period[(10^(-4))J],{J,1,2},AccuracyGoal->6,PrecisionGoal->6]
]

But:

ϕDefaultParametric=
ParametricNDSolveValue[
{ϕ''[t]==rhs/.defaultParameters,ϕ[0]==ϕInit,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,Sow[t];"StopIntegration"]},
ϕ,{t,0,100},{J,ϕInit}
]

Clear[periodParametric]
periodParametric[J_?NumericQ,ϕInit_?NumericQ,ϕ_ParametricFunction]:=Reap[ϕ[J,ϕInit]][[2,1,1]]

The following command works the first time only. After that, ParametricNDSolve uses caching, and the Sow is not executed again.

periodParametric[(10^(-4)),148°,ϕDefaultParametric]

periodParametric[(10^(-4)),148°,ϕDefaultParametric]

On 18.10.2024, Oliver Ruebenkoenig suggested to switch off caching. Then everything works:

ϕDefaultParametricNoCache=
ParametricNDSolveValue[
{ϕ''[t]==rhs/.defaultParameters,ϕ[0]==ϕInit,ϕ'[0]==0,
WhenEvent[ϕ'[t]<0,Sow[t];"StopIntegration"]},
ϕ,{t,0,100},{J,ϕInit},
Method->{"ParametricCaching"->None} (* Oliver Ruebenkoenig, 18.10.2024 *)
]

periodParametric[(10^(-4)),148°,ϕDefaultParametricNoCache]

Plot[periodParametric[(10^(-4)),ψ °,ϕDefaultParametric],{ψ,0.1,175},PlotStyle->Thick,PlotTheme->"Detailed",PlotLegends->None,BaseStyle->24,FrameLabel->{"initial angle [°]","period [s]"},ImageSize->800,Exclusions->None]

Plot[periodParametric[(10^(-4)),ψ °,ϕDefaultParametricNoCache],{ψ,0.1,175},PlotStyle->Thick,PlotTheme->"Detailed",PlotLegends->None,BaseStyle->24,FrameLabel->{"initial angle [°]","period [s]"},ImageSize->800,Exclusions->None]

Block[{T1=2,J1},FindRoot[period[(10^(-4))J1,90°,defaultParameters]==T1,{J1,1,2},AccuracyGoal->6,PrecisionGoal->6]]

Block[{T1=2,J1},FindRoot[periodParametric[(10^(-4))J1,90°,ϕDefaultParametricNoCache]==T1,{J1,1,2},AccuracyGoal->6,PrecisionGoal->6]]

Again, the previous command works only once with ϕDefaultParametric.

RepeatedTiming @ Block[{T1=2,J1},FindRoot[period[(10^(-4))J1,90°,defaultParameters]==T1,{J1,1,2},AccuracyGoal->6,PrecisionGoal->6]]

RepeatedTiming @ Block[{T1=2,J1},FindRoot[periodParametric[(10^(-4))J1,90°,ϕDefaultParametricNoCache]==T1,{J1,1,2},AccuracyGoal->6,PrecisionGoal->6]]

Appendix 4: DFixedPoints and DStabilityConditions

DFixedPoints[ϕ''[t]==-Sin[ϕ[t]],ϕ,t]

DStabilityConditions[ϕ''[t]==-Sin[ϕ[t]],ϕ,t]

DStabilityConditions[ϕ''[t]==rhs/.Rationalize[defaultParameters]/.J->(10^(-4)),ϕ,t]

DStabilityConditions[ϕ''[t]==-(ω^2)Sin[ϕ[t]],ϕ,t,Assumptions->ω>0]

DStabilityConditions[ϕ''[t]==rhs/.Rationalize[defaultParameters],ϕ,t,Assumptions->J>0]

AbsoluteTiming @ DStabilityConditions[ϕ''[t]==rhs,ϕ,t,Assumptions->m>0∧\[RawEscape]l>0∧r>0∧g>0∧J>0]

AbsoluteTiming @ DStabilityConditions[ϕ''[t]==rhs/.g->(981/100),ϕ,t,Assumptions->m>0∧\[RawEscape]l>2r>0∧J>0]

AbsoluteTiming @ Assuming[
m>0∧\[RawEscape]l>2r>0∧J>0,
DStabilityConditions[ϕ''[t]==rhs/.g->(981/100),ϕ,t]
]

The second end.

Preliminary material

SetOptions[EvaluationNotebook[],ScreenStyleEnvironment -> "Slideshow Working"]