Laboratorio di Fisica Computazionale
Capitolo 5. Problemi di Meccanica

5.1 Moto verticale nel campo della gravità

Senza resistenza dell'aria
Tutti i corpi si muovono con la stessa accelerazione costante
 g = 9.8 (m/(s^2))
F = m g  (forza peso di un corpo di massa m)
La legge oraria può essere ottenuta integrando l'equazione differenziale del II ordine

x''[t] = -g

oppure il sistema di due equazioni differenziali del I ordine

x'[t] = v[t]
v'[t] = -g

Clear[t,x,v,g,x0,v0];
x[t_]=x[t]/.DSolve[{x''[t] ==-g,x[0]==x0,x'[0]==v0},x[t],t][[1]]//Expand
v[t_]=x'[t]
t0=t/.Solve[x[t]==0,t][[2]]

Clear[t,x,v]
{x[t_],v[t_]}={x[t],v[t]}/.DSolve[{x'[t] == v[t],v'[t]==-g,x[0]==x0,v[0] == v0},{x[t],v[t]},t][[1]]//ExpandAll

val={g->9.8,b->1,x0->0,v0->98};
tmax=If[Im[t0/.val]==0,t0/.val,10];
Plot[x[t]/.val,{t,0,tmax},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t]/.val,{t,0,tmax},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

Con resistenza dell'aria
Introducendo la resistenza dell'aria
-m g - c v = m x''
x'[t] = v[t]
v'[t] = -g -b v[t]
b = c/m

Clear["Global`*"];
x[t_]=x[t]/.DSolve[{x''[t] +b  x'[t]==-g,x[0]==x0,x'[0]==v0},x[t],t][[1]]//Simplify
v[t_]=x'[t]
v[t]/.(E^(-b t))->0

val={g->9.8,b->1,x0->0,v0->98}
t0=t/.FindRoot[(x[t]/.val)==0,{t,(10/b)/.val}]
Plot[x[t]/.val,{t,0,t0},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"] 
Plot[{v[t]/.val},{t,0,t0},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

Con attrito radente
Introducendo la forza d'attrito radente 
(costante ed opposta al moto)

-mg  - c v/|v| = m x''

x'[t] = v[t]
v'[t] = -g -k v[t]/|v[t]| = - g -k sign(v)
k = c/m

Clear["Global`*"];
x1[t_]=x[t]/.DSolve[{x''[t]==-g -k Sign[v0] ,x[0]==x0,x'[0] == v0},x[t],t][[1]]//Simplify
tl=t/.Solve[x1'[t]==0,t][[1]];
x2[t_]=x[t]/.DSolve[{x''[t]==-g +k Sign[v0] ,x[0]==x1[tl],x'[0] == 0},x[t],t][[1]]//Simplify
x[t_]:=x1[t]/;t≤tl
x[t_]:=x2[t-tl]/;t>tl
v[t_]:=x1'[t]/;t≤tl
v[t_]:=x2'[t-tl]/;t>tl

val={g->9.8,k->5,x0->0,v0->100};
tl=tl/.val
t0=tl+t/.FindRoot[(x2[t]/.val)==0,{t,0.2`}];
Plot[x[t]/.val,{t,0,t0},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[{v[t]/.val},{t,0,t0},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

Clear[x,v];
{g,k,x0,v0}={9.8,5,0,98};
{x[t_],v[t_]}={x[t],v[t]}/.NDSolve[{x'[t] == v[t],v'[t]==-g -k Sign[v[t]],x[0]==x0,v[0] == v0},{x[t],v[t]},{t,0,t0}][[1]]

Plot[x[t],{t,0,t0},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t],{t,0,t0},PlotRange->All,AxesLabel->{"t (s)","v (m)"},PlotLabel->"Grafico velocità-tempo"]

5.2 Scivolamento di una catena appoggiata su un tavolo con attrito

Una catena di massa m e lunghezza a è appoggiata in parte su un tavolo orizzontale ed in parte pende verticalmete. Esaminare il moto nel campo della gravità
        a) senza attrito
        b) con attrito

Soluzione
Asse x verticale orientato verso il basso con origine al livello del tavolo.
x 		= tratto di cui pende la catena
a - x 	= tratto appoggiato sul tavolo

Forze agenti senza attrito
	F = mgx/a 	(0 ≤ x ≤ a)
		mg		(x > a)

Clear["Global`*"];
f[x_] = Which[x<a,g x/a,x≥a,g]
tmax = 4
{g,a,x0,v0}={9.8,10,2,0};
{x[t_],v[t_]}={x[t],v[t]}/.NDSolve[{x'[t] == v[t],v'[t]==f[x[t]],x[0]==x0,v[0] == v0},{x[t],v[t]},{t,0,tmax}][[1]]
Plot[x[t],{t,0,tmax},PlotRange->{0,40},AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t],{t,0,tmax},PlotRange->{0,30},AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]
Plot[v'[t],{t,0,tmax},PlotRange->{{0,tmax},{0,15}},AxesLabel->{"t (s)","a (m/(s^2))"},PlotLabel->"Grafico accelerazione-tempo"]

Forze agenti  con attrito

	F = 0 						(0 ≤ x < μa/(1+μ))
	F = mgx/a  - μmg(a-x)/a           (μa/(1+μ) ≤ x ≤ a)
	   =mg						(x > a)

Clear["Global`*"];
f[x_] =(g/a) Which[x<μ a/(1+μ),0,x<a,x - μ (a-x),x>a,a]
{g,a,μ,x0,v0}={9.8,10,0.2,2,0};
tmax=4
{x[t_],v[t_]}={x[t],v[t]}/.NDSolve[{x'[t] == v[t],v'[t]==f[x[t]],x[0]==x0,v[0] == v0},{x[t],v[t]},{t,0,tmax}][[1]]
Plot[x[t],{t,0,tmax},PlotRange->{0,50},AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t],{t,0,tmax},PlotRange->{0,30},AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]
Plot[v'[t],{t,0,tmax},PlotRange->{{0,tmax},{0,15}},AxesLabel->{"t (s)","a (m/(s^2))"},PlotLabel->"Grafico accelerazione-tempo"]

5.3 Moto verticale di un razzo

Impostazione del problema
Il motore di un razzo emette gas con velocità relativa vrel. Detta m la massa del razzo al tempo t, -dm la massa del gas emessa nell'intervallo di tempo tra t e t+dt, v la sua velocità e dv l'incremento di velocità subito nello stesso intervallo, la variazione di quantità di moto è data da

dp=(m-dm)(v+dv)+dm(v+vrel) - mv

dove le velocità sono prese con il loro segno. Trascurando i termini di ordine superiore al primo e dividendo per dt si ottiene

(dp/dt)= m(dv/dt)+vrel(dm/dt)

Dalla II legge della dinamica e considerando soltanto la forza peso

mx''[t] = -mg - m'[t] vrel
x''[t] = -g - vrel m'[t]/m[t]
deve essere specificata la funzione m[t]

Unità di misura
x in km
m in ton
g = 0.0098 km/(s^2)
mr = 6 ton (massa del razzo)
mc = 15 ton (massa del carburante)
rc = 0.19 ton/s (carburante bruciato al secondo)
vrel = 2.8 km/s (velocità di espulsione del carburante relativa al razzo)

Soluzione analitica (gravità costante)

Clear["Global`*"];
m[t_]=m0 - r*t;
sol=DSolve[x''[t] == - g - vrel m'[t]/m[t],x[t],t]//FullSimplify

Soluzione numerica (gravità costante)

Clear["Global`*"];
km=1;
ton=1;
s=1;
m=(km/(10^3));
kg=(ton/(10^3));
g=((9.8` m)/(s^2));
vrel=((2800 m)/s);
mr=6 ton;
mc=15 ton;
rc=((190 kg)/s);
tc=(mc/rc);
tmin=0 s;
tmax=150 s;
scal=(1/(10^3));
massa[t_]=If[t<tc,mr+mc-rc t,mr];
Plot[massa[t],{t,0,tmax},PlotRange->{0,mc+mr},AxesLabel->{"t (s)","massa (ton)"},PlotLabel->"Massa razzo+carburante"]
a[t_]=-g-((vrel massa'[t])/(massa[t]));
Plot[a[t],{t,tmin,tmax},AxesLabel->{"t (s)","a (km/\!\(\*SuperscriptBox[\(s\), \(2\)]\))"},PlotLabel->"accelerazione"]

Clear[sol,x]
tm=900 s;
sol=NDSolve[{x''[t]==-g-((vrel massa'[t])/(massa[t])),x[0]==0,x'[0]==0},x[t],{t,0,tm}];
x[t_]=x[t]/.sol[[1]];
tmax=t/.FindRoot[x[t]==0,{t,tm}];
Plot[x[t],{t,0,tmax},AxesLabel->{"t (s)","x (km)"},PlotLabel->"Grafico spazio-tempo"]
Plot[x'[t],{t,0,tmax},AxesLabel->{"t (s)","v (km/s)"},PlotLabel->"Grafico velocità-tempo"]

Soluzione numerica (gravità variabile)

Clear["Global`*"];
km=1;
ton=1;
s=1;
m=(km/(10^3));
kg=(ton/(10^3));
g=((9.8` m)/(s^2));
vrel=((2800 m)/s);
mr=6 ton;
mc=15 ton;
rc=((190 kg)/s);
tc=(mc/rc);
tm=1000 s;
massa[t_]=If[t<tc,mr+mc-rc t,mr];
RT=6400 km;
grav[x_]=If[x>0,(g/(((1+(x/RT)))^2)),g (1+(x/RT))];
sol=NDSolve[{x''[t]==-grav[x[t]]-((vrel massa'[t])/(massa[t])),x[0]==0,x'[0]==0},x[t],{t,0,tm}];
x[t_]=x[t]/.sol[[1]];
tmax=t/.FindRoot[x[t]==0,{t,tm}];
Plot[x[t],{t,0,tmax},AxesLabel->{"t (s)","x (km)"},PlotLabel->"Grafico spazio-tempo"]
Plot[x'[t],{t,0,tmax},AxesLabel->{"t (s)","v (km/s)"},PlotLabel->"Grafico velocità-tempo"]
fact=((10^3)/(3600^2));
a[t_]=fact (-grav[x[t]]-((vrel massa'[t])/(massa[t])));
Plot[a[t],{t,0,tmax},AxesLabel->{"t (s)","a (km/\!\(\*SuperscriptBox[\(s\), \(2\)]\))"},PlotLabel->"Grafico accelerazione-tempo"]
vfuga=(Sqrt[2 grav[x[tc]] (RT+x[tc])]);
Print["Velocità di fuga = ",vfuga];

5.4 Moti oscillatori

Oscillazioni armoniche
						F = -K x 			(ω^2)=(k/m)	
						
						 x''[t] +(ω^2) x[t]=0

Clear[sol,ω,x,v,x1,x2,x0,v0,rul1,rul2]
x1[t_]=x[t]/.DSolve[{x''[t]+(ω^2)x[t]==0},x[t],t][[1]]
x2[t_]=x1[t]/.{C[2]->A Cos[θ],C[1]->A Sin[θ]}//FullSimplify
rul1[xo_,vo_]=MapAt[Together,Solve[{x1[0]==xo,x1'[0]==vo},{C[1],C[2]}][[1]],{{1,2},{2,2}}]//Simplify
rul2[xo_,vo_]=MapAt[Together,Solve[{x2[0]==xo,x2'[0]==vo},{A,θ}][[1]],{{1,2},{2,2}}]//Simplify
x[t_]=x1[t]/.rul1[x0,v0]
v[t_]=x1'[t]/.rul1[x0,v0];

{m,ω,x0,v0}={2.5,0.314,10,0};
Plot[x[t],{t,0,40},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t],{t,0,40},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

k=m (ω^2);
T[t_]=(1/2)m ((v[t])^2);
V[t_]=(1/2)k ((x[t])^2);
Et[t_]=T[t]+V[t];

Plot[T[t],{t,0,40},PlotRange->{0,15},AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico Energia potenziale-tempo"]
Plot[V[t],{t,0,40},PlotRange->{0,15},AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico Energia cinetica-tempo"]
Plot[{T[t],V[t]/.rul1[x0,v0],Et[t]/.rul1[x0,v0]},{t,0,40},PlotRange->{0,15},AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico Energia totale-tempo"]

Oscillazioni smorzate
						F = -K x -b v			(ω^2)=(k/m)         δ=(b/(2 m))
						
						 x''[t] + 2 δ x'[t] + (ω^2) x[t]=0

Clear[sol,ω,Ω,δ,f,x,v,x1,x2,x3,x0,v0,rul1,rul2,rul3]
x1[t_,ω_,δ_]=x[t]/.DSolve[{x''[t]+2δ x'[t]+(ω^2)x[t]==0},x[t],t][[1]]
x2[t_,ω_,δ_]=ExpToTrig[x1[t,ω,δ]/.{(Sqrt[(δ^2)-(ω^2)])->I(Sqrt[(ω^2)-(δ^2)]),C[1]->A (E^(-I φ)),C[2]->A (E^(I φ))}]//FullSimplify
x3[t_,δ_]=x[t]/.DSolve[{x''[t]+2δ x'[t]+(δ^2)x[t]==0},x[t],t][[1]]//Simplify
rul1[ω_,δ_,x0_,v0_]=MapAt[Together,Solve[{x1[0,ω,δ]==x0,(D[x1[t,ω,δ],t]/.t->0)==v0},{C[1],C[2]}][[1]],{{1,2},{2,2}}]//Simplify
rul2[ω_,δ_,x0_,v0_]=MapAt[Together,Solve[{x2[0,ω,δ]==x0,(D[x2[t,ω,δ],t]/.t->0)==v0},{A,φ}][[4]],{{1,2},{2,2}}]//Simplify//PowerExpand
rul3[δ_,x0_,v0_]=MapAt[Together,Solve[{x3[0,δ]==x0,(D[x3[t,δ],t]/.t->0)==v0},{C[1],C[2]}][[1]],{{1,2},{2,2}}]//Simplify
x[t_,ω_,δ_,x0_,v0_]:=
Which[δ>ω,x1[t,ω,δ]/.rul1[ω,δ,x0,v0],
            δ<ω,x2[t,ω,δ]/.rul2[ω,δ,x0,v0],
           δ==ω,x3[t,δ]/.rul3[δ,x0,v0]
]
v[t_,ω_,δ_,x0_,v0_]:=
Which[δ>ω,((∂_x)x1[x,ω,δ]/.x->t)/.rul1[ω,δ,x0,v0],
            δ<ω,((∂_x)x2[x,ω,δ]/.x->t)/.rul2[ω,δ,x0,v0],
            δ==ω, ((∂_x)x3[x,δ]/.x->t)/.rul3[δ,x0,v0]
]

{x1[t,δ,δ],x2[t,δ,δ],x3[t,δ]}//Simplify

{Limit[rul1[ω,δ,x0,v0],ω->δ],Limit[rul2[ω,δ,x0,v0],ω->δ]}

{Limit[x1[t,ω,δ]/.rul1[ω,δ,x0,v0],ω->δ],Limit[x2[t,ω,δ]/.rul2[ω,δ,x0,v0],ω->δ],x3[t,δ]/.rul3[δ,x0,v0]}//Simplify//PowerExpand

Se δ >ω  la soluzione dell'omogenea decresce in modo esponenziale:

{ω,δ,x0,v0}={0.2,0.5,10,0};
Plot[{x[t,ω,δ,x0,v0]},{t,0,300 δ},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t,ω,δ,x0,v0],{t,0,300δ},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

Se  δ < ω  la soluzione dell'omogenea associata diventa oscillante:

{ω,δ,x0,v0}={0.657,0.1,10,0};
Plot[{x0 (E^(-δ t)),-x0 (E^(-δ t)),x[t,ω,δ,x0,v0]},{t,0,300 δ},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t,ω,δ,x0,v0],{t,0,300δ},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

Se  δ = ω  si ha lo smorzamento critico:

{ω,δ,x0,v0}={0.2,0.2,10,0};
Plot[x[t,ω,δ,x0,v0],{t,0,300δ},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t,ω,δ,x0,v0],{t,0,300δ},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]

Oscillazioni forzate
						F = -K x - b x'[t] + (Subscript[F, 0])Cos[Ω t]	(ω^2) =(k/m)   δ =(b/(2m))   f =((Subscript[F, 0])/m)
						
						 x''[t] +2 δ x'[t]+(ω^2) x[t] =f Cos[Ω t]

Clear[sol,ω,Ω,δ,f,x,x0,x01,x02,x03,xas,x1,x2,x3,a,b,A,B,F,φ,θ,xo,vo,xx,rul1,rul2]
sol=DSolve[x''[t]+2δ x'[t]+(ω^2)x[t]==f Cos[Ω t],x[t],t]

x01[t_]=x[t]/.sol[[1]]/.f->0
xas[t_]=x[t]/.sol[[1]]/.{C[1]->0,C[2]->0}//Simplify

Se δ < ω  la soluzione dell'omogenea associata diventa oscillante:

x02[t_]=ExpToTrig[x01[t]/.{(Sqrt[(δ^2)-(ω^2)])->I(Sqrt[(ω^2)-(δ^2)]),C[1]->A Sin[θ]+I A Cos[θ],C[2]->A Sin[θ]-I A Cos[θ]}]//FullSimplify

x1[t_]=x01[t]+xas[t]
x2[t_]=x02[t]+xas[t]

Calcoliamo le regole che permettono di ricavare  A e θ  o C[1] e C[2] dalla condizioni iniziali  x(0)=x0 v(0)=v0

Clear[xo,vo,δ,ω,Ω,f,xx]
rul1=MapAt[Together,Solve[{x1[0]==xo,x1'[0]==vo},{C[1],C[2]}][[1]],{{1,2},{2,2}}]//Simplify
sol=Solve[{x2[0]==xo,x2'[0]==vo}/.{ Sin[θ]->(u/A),Cos[θ]->(v/A)},{u,v}][[1]]//Simplify;
rul2={A->Sqrt[((u^2)+(v^2))/.sol//Simplify],θ->ArcTan[Numerator[(v/u)],Denominator[(v/u)]]/.sol//Simplify}

Calcoliamo le soluzioni per δ=0 e  δ = ω

xx[t_]=xx [t]/.DSolve[{xx''[t]+(ω^2)xx[t]==Cos[Ω t],xx[0]==xo,xx'[0]==vo},xx[t],t][[1]]//Simplify;
x3[t_]=Limit[x1[t]/.rul1,ω->δ]//FullSimplify
x03[t_]=Limit[x01[t]/.rul1,ω->δ]//FullSimplify

x0[t_]:=Which[δ==ω,x03[t],δ>ω,x01[t]/.rul1,δ<ω,x02[t]/.rul2];
x[t_,ω_,δ_/;δ>ω]:=x1[t]/.rul1
x[t_,ω_,δ_/;δ<ω]:=x2[t]/.rul2
x[t_,ω_,δ_/;δ==ω]:=x3[t]
v[t_,ω_,δ_/;δ>ω]:=x1'[t]/.rul1
v[t_,ω_,δ_/;δ<ω]:=x2'[t]/.rul2
v[t_,ω_,δ_/;δ==ω]:=x3'[t]

La soluzione asintotica può essere scrtitta nella forma:

{ω,δ,xo,vo,f,Ω}={1,0.1,0,1,1,Ω=2 ω};
Plot[x[t,ω,δ],{t,0,500δ},PlotRange->All,AxesLabel->{"t (s)","x (m)"},PlotLabel->"Grafico spazio-tempo"]
Plot[v[t,ω,δ],{t,0,500δ},PlotRange->All,AxesLabel->{"t (s)","v (m/s)"},PlotLabel->"Grafico velocità-tempo"]
Clear[ω,δ,xo,vo,f,Ω]

(a Sin[Ω t]+b Cos[Ω t])/.{a->F Cos[φ],b->F Sin[φ]}//FullSimplify

Con a e b dati da:

a=Simplify[Coefficient[xas[t],Sin[Ω t]]]
b=Simplify[Coefficient[xas[t],Cos[Ω t]]]

Ricaviamo infine F e φ

F[f_,δ_,ω_,Ω_]=PowerExpand[Simplify[(Sqrt[(a^2)+(b^2)])]]
φ[δ_,ω_,Ω_]=ArcTan[Numerator[(b/a)],Denominator[(b/a)]]

Se poniamo ω=1 e f=1 
				Ω    viene misurato in unità di ω 
				δ    viene misurato in unità di ω 
				F    in unità (f/(ω^2))

ω=1;f=1;
Plot[Evaluate[Table[F[f,δ,ω,Ω],{δ,0.1`,0.6`,0.1`}]],{Ω,0,2.5`},PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"\!\(\*FractionBox[\(Ω\), \(ω\)]\)","\!\(\*FractionBox[SuperscriptBox[\(ω\), \(2\)], \(f\)]\)F"}]
Plot[Evaluate[Table[φ[δ,ω,Ω],{δ,0.1`,0.6`,0.1`}]],{Ω,0,2.5`},PlotRange->{0,π},Ticks->{Automatic,{0,(π/2),π}},AxesLabel->{"\!\(\*FractionBox[\(Ω\), \(ω\)]\)","φ"}]