Programma che genera i numeri primi

prime={2};
Do[
Do[If[(k/(prime[[i]]))∈Integers,Goto[qua],0],
{i,1,Length[prime]}];
prime=Append[prime,k];
Label[qua],
{k,3,200}]
prime

Table[Prime[n],{n,46}]
prime==%

Animazioni di traiettorie

Clear[xx,yy]
xx[t_]=Cos[10*t];
yy[t_]=Sin[8*t];

(* solo particella *)Animate[Graphics[{PointSize[.03],Point[{xx[t],yy[t]}]},Axes->True,PlotRange->{{-1.1,1.1},{-1.1,1.1}}],{t,0,5},AnimationRunning->False]

(* solo traiettoria *)Animate[ParametricPlot[{xx[t],yy[t]},{t,0,a},PlotRange->{{-1.1,1.1},{-1.1,1.1}}],{a,0,5},AnimationRunning->False]

(* particella con traiettoria  *)
Animate[Show[
ParametricPlot[{xx[t],yy[t]},{t,0,a},PlotRange->{{-1.1,1.1},{-1.1,1.1}}],Graphics[{PointSize[.03],Point[{xx[a],yy[a]}]},Axes->True,PlotRange->{{-1.1,1.1},{-1.1,1.1}}]],{a,0,5},AnimationRunning->False]

Scivolamento di una catena appoggiata su un tavolo

Una catena di massa m e lunghezza a è appoggiata in parte su un tavolo orizzontale ed in parte pende verticalmente. Esaminare il moto nel campo della gravità

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

	F = mgx[t]/a 	(0 ≤ x[t] ≤ a)
		mg		(x[t] > a)
		
equazione del moto:

m*x’’[t]==F

Clear[f,x,t,v];
f[x_] = Which[x<a,g x/a,x≥a,g];
tmax = 3;
{g,a,x0,v0}={9.8,.8,.01,0};
sol=NDSolve[{x''[t] == f[x[t]],x[0]==x0,x'[0] == v0},x[t],{t,0,tmax}]

Clear[x]
x[t_]=x[t]/.sol[[1]]

Plot[x[t],{t,0,1.7}]

t1=t/.FindRoot[x[t]==a,{t,1}]

Animate[Graphics[{Red,Rectangle[{0,0},{1,1}],Black,Thickness[.01],Line[{{1-(a-x[t]),1},{1,1}}],Line[{{1,1},{1,1-x[t]}}]},Axes->True,PlotRange->{{0,2},{0,1.5}}],{t,0,t1},AnimationRunning->False,AnimationRepetitions->1]

Problema di n corpi

Le eq. del moto  sono:

(Subscript[(OverVector[f]), i])=(Subscript[m, i]) (∂_(t,t))(Subscript[(OverVector[x]), i])       ->      (Subscript[f, i,j])=(Subscript[m, i]) (∂_(t,t))(Subscript[x, i,j])    

(l'indice i numera i corpi, l'indice j numera le 3 componenti spaziali)

Remove[fvec,x,m,cost]
n=4; (* numero dei corpi *)
(* ESEMPIO: forze centrali tipo (1/(r^2)):   *)
fvec[i_]:=Table[
Sum[If[i==j,0,((cost*(Subscript[m, i])*(Subscript[m, j]))/((((Sqrt[(∑_(p=1)^3)((((Subscript[x, j,p])[t]-(Subscript[x, i,p])[t]))^2)])))^3)) *((Subscript[x, j,k])[t]-(Subscript[x, i,k])[t])],{j,1,n}],
{k,1,3}](* =(Subscript[(OverVector[f]), i])  *)
eqvec=Flatten[Table[fvec[i][[j]]-(Subscript[m, i])*(∂_(t,t))(Subscript[x, i,j])[t]==0,{i,1,n},{j,1,3}]];(* eq. del moto  *)

inivec=Join[
Table[(Subscript[x, i,j])[0]==0,{i,1,n},{j,1,3}],
Table[(Subscript[x, i,j])'[0]==0,{i,1,n},{j,1,3}]
]//MatrixForm

inivec=Flatten[(({{(Subscript[x, 1,1])[0]==13, (Subscript[x, 1,2])[0]==0, (Subscript[x, 1,3])[0]==0}, {(Subscript[x, 2,1])[0]==12, (Subscript[x, 2,2])[0]==0, (Subscript[x, 2,3])[0]==0}, {(Subscript[x, 3,1])[0]==0, (Subscript[x, 3,2])[0]==0, (Subscript[x, 3,3])[0]==0}, {(Subscript[x, 4,1])[0]==3, (Subscript[x, 4,2])[0]==3, (Subscript[x, 4,3])[0]==2}, {(Subscript[x, 1,1])'[0]==0, (Subscript[x, 1,2])'[0]==.9, (Subscript[x, 1,3])'[0]==.3}, {(Subscript[x, 2,1])'[0]==0, (Subscript[x, 2,2])'[0]==1.1, (Subscript[x, 2,3])'[0]==0}, {(Subscript[x, 3,1])'[0]==0, (Subscript[x, 3,2])'[0]==0, (Subscript[x, 3,3])'[0]==0}, {(Subscript[x, 4,1])'[0]==-1.7, (Subscript[x, 4,2])'[0]==0, (Subscript[x, 4,3])'[0]==0}}))];(* condizioni iniziali  *)
cost=.01;
mvec={1,10,1.15*(10^3),.05};(* masse *)
Do[(Subscript[m, i])=mvec[[i]],{i,1,n}];
all=Join[eqvec,inivec];
incvec=Flatten[Table[(Subscript[x, i,j])[t],{i,1,n},{j,1,3}]];
tmax=200;
sol=Flatten[NDSolve[all,incvec,{t,0,tmax}]] (* soluzione *)
Do[
(Subscript[x, i,j])[t_]=(Subscript[x, i,j])[t]/.sol,
{i,1,n},{j,1,3}]

range=20;
Animate[Graphics[{
PointSize[.02],Point[{(Subscript[x, 1,1])[t],(Subscript[x, 1,2])[t]}],
Blue,PointSize[.03],Point[{(Subscript[x, 2,1])[t],(Subscript[x, 2,2])[t]}],
Red,PointSize[.04],Point[{(Subscript[x, 3,1])[t],(Subscript[x, 3,2])[t]}],
Blue,PointSize[.01],Point[{(Subscript[x, 4,1])[t],(Subscript[x, 4,2])[t]}]
},Axes->True,PlotRange->{{-range,range},{-range,range}}],
{t,0,tmax},AnimationRunning->False]

range=20;
Animate[Graphics3D[{
Red,Sphere[{(Subscript[x, 1,1])[t],(Subscript[x, 1,2])[t],(Subscript[x, 1,3])[t]},.3],
Blue,Sphere[{(Subscript[x, 2,1])[t],(Subscript[x, 2,2])[t],(Subscript[x, 2,3])[t]},.8],
Yellow,Sphere[{(Subscript[x, 3,1])[t],(Subscript[x, 3,2])[t],(Subscript[x, 3,3])[t]},3],
White,Sphere[{(Subscript[x, 4,1])[t],(Subscript[x, 4,2])[t],(Subscript[x, 4,3])[t]},.3]
},Boxed->False,Background->Black,PlotRange->{{-range,range},{-range,range},{-range,range}}],{t,0,tmax},AnimationRunning->False]

(* Conservazione della quantita' di moto totale *)
Plot[(∑_(i=1)^n)(Subscript[m, i])*{(Subscript[x, i,1])'[t],(Subscript[x, i,2])'[t],(Subscript[x, i,3])'[t]},{t,0,tmax},PlotRange->All]

(* Conservazione energia totale *)
Clear[kin,pot,ener]
kin[t_]=(∑_(i=1)^n)(1/2)*(Subscript[m, i])*((((Subscript[x, i,1])'[t])^2)+(((Subscript[x, i,2])'[t])^2)+(((Subscript[x, i,3])'[t])^2));
pot[t_]=(1/2)*Sum[If[i==j,0,-((cost*(Subscript[m, i])*(Subscript[m, j]))/(Sqrt[(∑_(k=1)^3)((((Subscript[x, j,k])[t]-(Subscript[x, i,k])[t]))^2)]))],{i,1,n},{j,1,n}];
ener[t_]=kin[t]+pot[t];
Plot[{kin[t],pot[t],ener[t]},{t,0,tmax}]

Corpo rigido con un punto fisso
Giroscopio

Studiamo il moto di un insieme di N particelle che hanno distanze fisse tra loro e dall’origine di un sistema inerziale x’ y’ z’ (corpo rigido con un punto fisso). L’insieme ha tre gradi di liberta’: i due angoli sferici θ(t) e φ(t) del versore (OverVector[k]) di una terna ortonormale {(OverVector[i]),(OverVector[j]),(OverVector[k])} solidale col corpo costruita come in figura, e l'angolo di rotazione α(t) intorno a (OverVector[k]).
Esprimiamo il vettore posizione (Subscript[(OverVector[x]), i]) della particella i in termini di θ[t],φ[t],α[t]:

(OverVector[k])={Sin[θ[t]]*Cos[φ[t]],Sin[θ[t]]*Sin[φ[t]],Cos[θ[t]]};
(OverVector[u])=((OverVector[k])/(Sin[θ[t]]))×{0,0,1};
(OverVector[n])=(OverVector[u])×(OverVector[k]);
(OverVector[i])=Cos[α[t]]*(OverVector[n])+Sin[α[t]]*(OverVector[u]);
(OverVector[j])=-Sin[α[t]]*(OverVector[n])+Cos[α[t]]*(OverVector[u]);
(Subscript[(OverVector[x]), i])=(Subscript[x, i])*(OverVector[i])+(Subscript[y, i])*(OverVector[j])+(Subscript[z, i])*(OverVector[k])  ==  (Subscript[(OverVector[x]), i])(θ[t],φ[t],α[t])

dove {(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])} sono le coordinate della particella i nel sistema {(OverVector[i]), (OverVector[j]),(OverVector[ k])} solidale col corpo.
(Subscript[(OverVector[x]), i]) e' una forma lineare nelle {(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])}: dunque il momento angolare (Subscript[(OverVector[j]), i]) e' una forma bilineare:

(Subscript[(OverVector[j]), i])=(Subscript[m, i])*(Subscript[(OverVector[x]), i])×(∂_t) (Subscript[(OverVector[x]), i])==Φ.(Subscript[(OverVector[μ]), i])  

dove:

(Subscript[(OverVector[μ]), i])= (Subscript[m, i])*{(Subscript[x, i]^2),(Subscript[x, i]) (Subscript[y, i]),(Subscript[y, i]^2),(Subscript[x, i]) (Subscript[z, i]),(Subscript[y, i]) (Subscript[z, i]),(Subscript[z, i]^2)} (momenti d’inerzia);
Φ = matrice 3×6 di funzioni delle  θ[t], φ[t], α[t] indipendente dal punto i.

Le 3 eq. del moto nelle 3 variabili  θ[t],φ[t],α[t] sono:

(Subscript[(OverVector[τ]), est])=(∂_t) (OverVector[J])       ->     (∑_(i=1)^N)(Subscript[(OverVector[x]), i])×(Subscript[(OverVector[f]), i, est]) =(∂_t) (∑_(i=1)^N)(Subscript[(OverVector[j]), i])      ->       (∑_(i=1)^N)(Subscript[(OverVector[x]), i])×(Subscript[(OverVector[f]), i,est]) =  ((∂_t) Φ).( (∑_(i=1)^N) (Subscript[(OverVector[μ]), i]) )

dove  (Subscript[(OverVector[f]), i est]) e' la forza esterna sulla particella i; nel caso di corpi continui  (∑_(i=1)^N) (Subscript[m, i])  ->  ∫dm .

-----------------------------------------------------------------

Esercizio: verificare le relazioni precedenti.

Esercizio: studiare il moto di un giroscopio con: sola forza di gravita’; gravita’ con attrito dell’aria; gravita’ con giroscopio accelerato da un momento motore.

------------------------------------------------------------------

(*  vettore posizione (Subscript[(OverVector[x]), i]) della particella i in termini di θ[t],φ[t],α[t]   *)
Remove[x,y,z,k,u,n,i,j,θ,φ,α,m,η,μ];
(OverVector[k])={Sin[θ[t]]*Cos[φ[t]],Sin[θ[t]]*Sin[φ[t]],Cos[θ[t]]};
(OverVector[u])=((OverVector[k])/(Sin[θ[t]]))\[Cross]{0,0,1};
(OverVector[n])=(OverVector[u])\[Cross](OverVector[k])//Simplify;
(OverVector[i])=Cos[α[t]]*(OverVector[n])+Sin[α[t]]*(OverVector[u]);
(OverVector[j])=-Sin[α[t]]*(OverVector[n])+Cos[α[t]]*(OverVector[u]);
(Subscript[(OverVector[x]), i_])=(Subscript[x, i])*(OverVector[i])+(Subscript[y, i])*(OverVector[j])+(Subscript[z, i])*(OverVector[k])//Simplify;
MatrixForm[%]

(*  Momento angolare della particella i: (Subscript[(OverVector[j]), i])(t) = Φ(t).(Subscript[(OverVector[μ]), i])  *)
(Subscript[(OverVector[η]), i_])={(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])};
(Subscript[(OverVector[μ]), i_])=(Subscript[m, i])*Flatten[Table[Table[(Subscript[(OverVector[η]), i])[[ii]]*(Subscript[(OverVector[η]), i])[[jj]],{ii,1,jj}],{jj,1,3}]];
(Subscript[(OverVector[j]), i_])=(Subscript[m, i])*(Subscript[(OverVector[x]), i])\[Cross](∂_t) (Subscript[(OverVector[x]), i]);
Φ=Table[Coefficient[(Subscript[(OverVector[j]), i])[[ii]],(Subscript[(OverVector[μ]), i])[[jj]]],{ii,1,3},{jj,1,6}]//Simplify;
((Subscript[(OverVector[j]), i])-Φ.(Subscript[(OverVector[μ]), i]))//Simplify
(* (Subscript[(OverVector[j]), i]) e' una forma lineare nei momenti d'inerzia (Subscript[(OverVector[μ]), i])={(Subscript[m, i]) (Subscript[x, i]^2),(Subscript[m, i]) (Subscript[x, i]) (Subscript[y, i]),(Subscript[m, i]) (Subscript[y, i]^2),(Subscript[m, i]) (Subscript[x, i]) (Subscript[z, i]),(Subscript[m, i]) (Subscript[y, i]) (Subscript[z, i]),(Subscript[m, i]) (Subscript[z, i]^2)} *)

(* Soluzione e animazione delle equazioni del moto *)
(* Giroscopio *)
(({{(Subscript[m, 1]), (Subscript[x, 1]), (Subscript[y, 1]), (Subscript[z, 1])}, {(Subscript[m, 2]), (Subscript[x, 2]), (Subscript[y, 2]), (Subscript[z, 2])}, {(Subscript[m, 3]), (Subscript[x, 3]), (Subscript[y, 3]), (Subscript[z, 3])}, {(Subscript[m, 4]), (Subscript[x, 4]), (Subscript[y, 4]), (Subscript[z, 4])}}))=(({{10, .3, 0, 1}, {10, 0, .3, 1}, {10, -.3, 0, 1}, {10, 0, -.3, 1}})); ;(* coordinate delle particelle nel sistema solidale *)
g=9.8;
Clear[θ,φ,α];
Remove[f];
(Subscript[(OverVector[f]), i_, est])={0,0,-(Subscript[m, i])*g};
num=4;
vvv=((∑_(i=1)^num)(Subscript[(OverVector[x]), i])\[Cross](Subscript[(OverVector[f]), i, est]) )-((∂_t) Φ).( (∑_(i=1)^num) (Subscript[(OverVector[μ]), i]) );(* equazioni del moto *)
eqvec=Table[vvv[[i]]==0,{i,1,3}];
all=Join[eqvec,{θ[0]==(π/4),φ[0](* ≠0 *)==(π/2),α[0](* =0 *)==0,θ'[0]==0,φ'[0]==0,α'[0]==120}]; (* vettore delle equazioni piu' le condizioni iniziali *)
tmax=10;
sol=Flatten[NDSolve[all,{θ[t],φ[t],α[t]},{t,0,tmax},AccuracyGoal->6,PrecisionGoal->6,MaxSteps->50000]];
θ[t_]=θ[t]/.sol
φ[t_]=φ[t]/.sol
α[t_]=α[t]/.sol

range=2.1;
Animate[Graphics3D[{
Sphere[(Subscript[(OverVector[x]), 1])/.t->tt,.2],
Sphere[(Subscript[(OverVector[x]), 2])/.t->tt,.2],
Sphere[(Subscript[(OverVector[x]), 3])/.t->tt,.2],
Sphere[(Subscript[(OverVector[x]), 4])/.t->tt,.2],
Thickness[.004],
Line[{{0,0,0},2*(OverVector[k])/.t->tt}]
},Axes->True,PlotRange->{{-range,range},{-range,range},{-range,range}}],{tt,0,tmax},AnimationRunning->False]

(* precessione e nutazione *)
range=2.1;
Show[ParametricPlot3D[2*(OverVector[k])/.t->tt,{tt,0,tmax},PlotStyle->Thickness[.003],PlotRange->{{-range,range},{-range,range},{-range,range}}],Graphics3D[{Arrow[{{0,0,0},2*(OverVector[k])/.t->0}],Opacity[.4],Sphere[{0,0,0},2]},Axes->True,PlotRange->{{-range,range},{-range,range},{-range,range}}]]

-------------------------------------------------------------

Esercizio: determinare il vettore velocita’ angolare in funzione di θ[t], φ[t], α[t].

La velocita’ di ogni punto i di un corpo rigido e’ data da:

(∂_t)(Subscript[(OverVector[x]), i])= (∂_t)(Subscript[(OverVector[x]), o]) +  (OverVector[ω]) × (Subscript[(OverVector[x]), i])

Nel caso di corpo rigido con un punto fisso ho (∂_t)(Subscript[(OverVector[x]), o])=0 ed ottengo (OverVector[ω]) risolvendo la precedente eq. su due punti generici (Subscript[(OverVector[x]), 1]) e (Subscript[(OverVector[x]), 2]):

Remove[ω]
(OverVector[ω])={(Subscript[ω, 1]),(Subscript[ω, 2]),(Subscript[ω, 3])};
Solve[{((∂_t) (Subscript[(OverVector[x]), 1]))==(OverVector[ω])\[Cross](Subscript[(OverVector[x]), 1]),((∂_t) (Subscript[(OverVector[x]), 2]))==(OverVector[ω])\[Cross](Subscript[(OverVector[x]), 2])},(OverVector[ω])]//Simplify;
MatrixForm[%[[1]]]

Esercizio: calcolare i momenti d’inerzia di un corpo rigido continuo (ellissoide, cubo...) con un punto fisso e studiarne il moto.

Esercizio: estendere il programma precedente al caso del moto di un corpo rigido libero, aggiungendo le altre tre equazioni relative al moto del centro di massa.

---------------------------------------------------------------