Foucaultovo nihalo

Enačbe gibanja Foucaultovega nihala so sistem diferencialnih enačb 

m D[ζ,{t,2}] == ( -m g {0,0,1} + T - 2m Cross[ω,D[ζ,t]])

z algebraičnim pogojem ‖ζ - l{0,0,1}‖== l.

Brezdimenzijski zapis

Ker bomo reševali numerično, je ugodno enačbe gibanja zapisati v brezdimenzijski obliki. Dimenzija časa je (Sqrt[(l/g)]) , enotska dolžina je l.

Krajevni vektor do točke nihala

ζ = {l ζ1[(Sqrt[(g/l)])t],l ζ2[(Sqrt[(g/l)])t],l ζ3[(Sqrt[(g/l)])t]}

Tu so ζ1, ζ2, ζ3 funkcije  brezdimenzijskega časa. Vektor kotne hitrosti zapisan v RKS je

ω = (Sqrt[(g/l)])ω0{-Sin[θ],0,Cos[θ]}

Tu je θ geografska širina, RKS pa je postavljen tako, da je os ζ3 v smeri navpičnice, ζ2 pa kaže proti vzhodu. 
Sila vrvice je dana z

vecT = ((T0[τ] m g)/l){-l ζ1[(Sqrt[(g/l)])t],-l ζ2[(Sqrt[(g/l)])t],l-l ζ3[(Sqrt[(g/l)])t]}

Enačba gibanja

eq=Simplify[(1/(m g))(m D[ζ,{t,2}]-( -m g {0,0,1} + vecT - 2m Cross[ω,D[ζ,t]]))]/.(Sqrt[(g/l)]) t-> τ

Enačba vezi

cnst= {-ζ1[τ],-ζ2[τ],1-ζ3[τ]}.{-ζ1[τ],-ζ2[τ],1-ζ3[τ]}== 1

Začetni pogoj, nihalo odklonemo v smer osi ζ1 in spustimo.

zp = {ζ1[0]== ζ0,ζ2[0]== 0,ζ3[0]== 1-(Sqrt[1-((ζ0)^2)]),ζ1'[0]== 0,ζ2'[0]== 0,ζ3'[0]== 0}

eqeq=Join[Thread[eq== 0],{cnst},zp]

Reševanje algebraično diferencialnega sistema

Ker rešujemo numerično potrebujemo konkretne podatke.

primer ={ ω0->(10^(-2)),ζ0-> 1/50,θ-> π/4}

tmax=200

Rešujemo algebraično diferencialno enačbo. Potrebno je sodelovanje uporabnika. Opcija Method->{“IndexReduction”->{Automatic,”ConstraintMethod”->None}} pripelje do rešitve.

slv=NDSolve[eqeq/.primer,{ζ1,ζ2,ζ3,T0},{τ,0,tmax},Method->{"IndexReduction"->{Automatic,"ConstraintMethod"->None}},MaxSteps->50000 ]

Preverimo, če je Mathematica pravilno rešila sistem.

Plot[Evaluate[eq/.primer/.slv],{τ,0,1},PlotRange->{Automatic,{-(10^(-4)),(10^(-4))}}]

Plot[Evaluate[((ζ1[τ])^2)+((ζ2[τ])^2)+(((1-ζ3[τ]))^2)/.slv],{τ,0,tmax},PlotRange->All]

Grafična predstavitev rešitve

ParametricPlot[Evaluate[{ζ1[τ],ζ2[τ]}/.slv],{τ,0,tmax},PlotRange->All]

slika1=Plot[Evaluate[{ζ2[t],ζ1[t]}/.slv[[1]]],{t,0,tmax},PlotRange->All]

Aproksimacija ζ3 = ζ3' = ζ3'' = 0

aprx={ζ3[τ]-> 0,(ζ3)'[τ]-> 0,(ζ3)''[τ]-> 0}

Eliminacija sile vrvice

slvT=Solve[(eq[[3]]/.aprx)== 0,T0[τ]]

eqa=eq[[1;;2]]/.slvT[[1]]/.aprx

zpa = {ζ1[0]== ζ0,ζ2[0]== 0,ζ1'[0]== 0,ζ2'[0]== 0}

Dobili smo nelinearen sistem diferencialnih enab

slva=NDSolve[Join[Thread[eqa== 0],zpa]/.primer,{ζ1,ζ2},{τ,0,tmax},MaxSteps->50000 ]

Plot[Evaluate[eqa/.primer/.slva],{τ,0,1},PlotRange->{Automatic,{-(10^(-4)),(10^(-4))}}]

ParametricPlot[Evaluate[{ζ1[τ],ζ2[τ]}/.slva],{τ,0,tmax},PlotRange->All]

slika2=Plot[Evaluate[{ζ2[t],ζ1[t]}/.slva[[1]]],{t,0,tmax},PlotRange->All]

Linearizacija, harmonična aproksimacija

Eliminacija nelinearnih členov

eqh={-2 ω0 Cos[θ] (ζ2)'[τ]+ζ1[τ] +(ζ1)''[τ],2 ω0 Cos[θ] (ζ1)'[τ]+ζ2[τ] +(ζ2)''[τ]}

zpa = {ζ1[0]== ζ0,ζ2[0]== 0,ζ1'[0]== 0,ζ2'[0]== 0}

slvh=NDSolve[Join[Thread[eqh== 0],zpa]/.primer,{ζ1,ζ2},{τ,0,tmax},MaxSteps->50000 ]

ParametricPlot[Evaluate[{ζ1[τ],ζ2[τ]}/.slvh],{τ,0,tmax},PlotRange->All]

slika3=Plot[Evaluate[{ζ2[t],ζ1[t]}/.slvh[[1]]],{t,0,tmax},PlotRange->All]

Primerjava

Show[{slika1,slika2,slika3}]

Vse tri metode dajo približno enak rezultat

Primerjava rešitev

Majhen odmik

primer ={ ω0->1/10,ζ0-> 1/50,θ-> π/4}

tmax=200

slv=NDSolve[eqeq/.primer,{ζ1,ζ2,ζ3,T0},{τ,0,tmax},Method->{"IndexReduction"->{Automatic,"ConstraintMethod"->None}},MaxSteps->50000 ]

slva=NDSolve[Join[Thread[eqa== 0],zpa]/.primer,{ζ1,ζ2},{τ,0,tmax},MaxSteps->50000 ]

slvh=NDSolve[Join[Thread[eqh== 0],zpa]/.primer,{ζ1,ζ2},{τ,0,tmax},MaxSteps->50000 ]

Plot[Evaluate[{ζ1[t]/.slv[[1]],ζ1[t]/.slva[[1]],ζ1[t]/.slvh[[1]]}],{t,0,tmax},PlotRange->All,PlotStyle-> {Red,Blue,Green}]

Velik odmik

primer ={ ω0->1/10,ζ0-> 1/5,θ-> π/4}

tmax=200

slv=NDSolve[eqeq/.primer,{ζ1,ζ2,ζ3,T0},{τ,0,tmax},Method->{"IndexReduction"->{Automatic,"ConstraintMethod"->None}},MaxSteps->50000 ]

slva=NDSolve[Join[Thread[eqa== 0],zpa]/.primer,{ζ1,ζ2},{τ,0,tmax},MaxSteps->50000 ]

slvh=NDSolve[Join[Thread[eqh== 0],zpa]/.primer,{ζ1,ζ2},{τ,0,tmax},MaxSteps->50000 ]

Plot[Evaluate[{ζ1[t]/.slv[[1]],ζ1[t]/.slva[[1]],ζ1[t]/.slvh[[1]]}],{t,0,tmax},PlotRange->All,PlotStyle-> {Red,Blue,Green}]

Zaključek

Za majhne pomike je linearizacija dobra aproksimacija