Problemi di Fisica con l'uso del PC

X Lezione

Procedure
Una procedura è una sequenza di espressioni separate da ";". Mathematica valuta le espressioni consecutivamente e restituisce il risultato dell'ultima espressione.

r  =  3;  s  =  Sqrt[4]  +  r;  r^2  +  s^2

Se non si vuole il risultato dell'ultima espressione bisogna farla seguire da un ";"

Print["What a wonderful world!"];
                                                     Plot[
        Sin[x],  {x,  0,  2 Pi}];
                                                     x  =  y;

Con una procedura possiamo costruire una funzione. In questo caso occorre usare delle parentesi "(.....)" per delimitare il corpo della funzione.

\(func[
        y_]\  := \ \ \((Print["What a wonderful world!"]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Plot[
          Sin[x], \ {x, \ 0, \ 
            2\ Pi}]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\  = \ y;)\)\)

Un effetto indesiderato è l'introduzione di una serie di variabili che non interessano

{r,  s,  x}

ClearAll[func,  r,  s,  x]

Simboli locali
Mathematica presenta un costrutto che consente di definire simboli che sono locali alla procedura e non vengono vengono visti dall'esterno. 
Module[{x,y,z=z0,....},procedura]   fa eseguire procedura, ma le variabili x,y,z,... sono solo locali e alla variabile z è assegnato inizialmente il valore z0.
Mathematica assegna alla variabile simbolo contenuta nel modulo il nome simbolo$n, dove n viene incrementato di uno per ogni chamata di Module in modo da distinguerle da altre variabili globali con lo stesso nome.

p  =  10
\(Module[{q, \ p}, \ Print[q]; \ q\  = \ p^2; \ 
      Print[q, \ "
"]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ p\  = \ 50; \ 
      Print[p, \ "   ", \ q]]\)
p
q
\(Module[{q, \ p}, \ Print[q]; \ q\  = \ p^2; \ 
      Print[q, \ "
"]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ p\  = \ 50; \ 
      Print[p, \ "   ", \ q]]\)

Notare che nell'inizializzazione delle variabili viene usata la variabile globale.

u  =  5;
Module[{u  =  120,  v  =  u^2  +  1},  Print[v]]
ClearAll[p,  u]

Funzioni condizionali
Una funzione condizionale restituisce il valore di una tra possibili espressioni alternative.
Le funzioni condizionali sono:
If[condizione,t]				restituisce t se condizione è vera, altrimenti restituisce Null
If[condizione,t,f]			restituisce t se condizione è vera,  f se condizione è falsa
If[condizione,t,f,u]			restituisce t se condizione è vera,  f se condizione è falsa, u se condizione 
					non è nè vera, nè falsa
Which[condizione-1,valore-1,		valuta in sequenza le condizion-i e restituisce il valore-i della prima che è vera.
condizione-2,valore-2,....]
Swhich[espressione,forma-1,valore-1,	valuta espressione e la confronta con le forme-i, appena trova uguaglianza  
forma-2,valore-2,....]			restituisce valore-i

Esempio: definizione della delta di Kronecker

kDel[n_Integer,  m_Integer]  :=  If[n  ==  m,  1,  0]
{kDel[1,  1],  kDel[1,  0]}

Esempio: definizione alternativa

newkDel[n_Integer,  m_Integer]  :=  1  /;  n  ==  m
newkDel[n_Integer,  m_Integer]  :=  0  /;  n  !=  m
{newkDel[1,  1],  newkDel[1,  0]}
ClearAll[kDel, newkDel]

Esempio: funzioni a gradino

\(v[x_]\  := \ 
      Which[x\  < \ \(-3\)\  || \ x\  >= \ 3, \ 
        0, 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((x\  >= \ \(-3\)\  && \ x\  < \ \(-1\))\) || \((x\  >= \ 1\  && \ 
              x\  < \ 3)\), \ \(-2\), 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\  >= \ \(-1\)\  && \ 
          x\  < \ 1, \ \(-1\)]\)
\(\(Plot[v[x], \ {x, \ \(-5\), \ 5}, \ 
        PlotStyle\  -> \ 
          Thickness[
            0.01], 
\ \ \ \ \ \ \ \ \ \ \ \ AxesLabel\  -> \ {" x", \ "v"}];\)\)

Definizione alternativa di funzione a gradino

newv[x_  /;  x  <  -3  ||  x  >=  3]  :=  0
newv[
        x_  /;  (x  >=  -3  &&  
                x  <  -1) || (x  >=  1  &&  
                x  <  3)]  :=  -2
newv[
        x_  /;  x  >=  -1  &&  x  <  1]  :=  -1

Notare che la prima definizione è differenziabile, la seconda no.

newv'[x]
newv'[1]
v'[x]
v'[1]
ClearAll[v,  newv]

Moto di un pianeta

(afelio e perielio misurati in milioni di Chilometri)
Si consideri il moto di un pianeta di massa m dovuto alla forza di attrazione da parte del Sole (massa M)
					F = G (( m M)/(r^2))
dove r è il modulo del vettore posizione del pianeta nel sistema di riferimento del Sole e G la costante di gravitazione universale. 
In forma vettoriale la II Legge della Dinamica si scrive: 
			\(TraditionalForm\`F\&_ = \(\(-G\) \(\(m\ M\)\/r\^3\) 
            r\&_ = \(m  a\&_ = m \(\( d\^2\) r\&_\)\/\(d\ t\^2\)\)\)\)
Proiettando questa equazione sui due assi di un sistema di riferimento cartesiano si ha
				\(\(\(d\^2\) x\)\/\(d\ t\^2\) + \(G\ M\ x\)\/\((x\^2 + y\^2)\)\^\(3/2\)\)=0
				\(\(\(d\^2\) y\)\/\(d\ t\^2\) + \(G\ M\ y\)\/\((x\^2 + y\^2)\)\^\(3/2\)\)=0
Se si usa un sistema di unità di misura in cui l'unità di lunghezza è l'Unità Astronomica, cioè la lunghezza del semiasse maggiore dell'orbita (ellittica) della Terra, e l'unità di tempo l'anno, si trova G M=4(π^2).

Disegnare la traiettoria di un pianeta che al tempo t=0 ha posizione e velocità date da
				x(0)=1, y(0)=0, x'(0)=-π, y'(0)=2π

\(
Clear[orbit]\)
orbit = 
      NDSolve[{x''[
                t] + (4 Pi^2) x[
                    t]/((x[t]^2 + y[t]^2)^(3/2)) == 
            0, y''[
                t] + (4 Pi^2) y[
                    t]/((x[t]^2 + y[t]^2)^(3/2)) == 0, 
          x[0] == 1, 
          y[0] == 0, x'[0] == -Pi, y'[0] == 
            2 Pi}, {x, y}, {t, 0, 1.6}]
ParametricPlot[Evaluate[{x[t], y[t]} /. orbit], {t, 0, 1.6}, 
        AspectRatio -> Automatic, 
        AxesLabel -> {x (AU), y (AU)}];

Ovviamente per disegnare un orbita occorre conoscere il periodo. Lo ricaviamo nel caso più semplice di orbite circolari. In questo caso la forza è tutta forza centripeta:

				F = G (( m M)/(r^2)) = m r (ω^2)
e, poiché ω = ((2  π)/T), 
				(T^2) = ((4  (π^2))/(G M)) (r^3)
Nel caso di orbite ellittiche questa formula diventa

				(T^2) = ((4  (π^2))/(G M)) (a^3)
dove a è il semiasse maggiore dell'ellissi.
Quindi per ricavare il periodo è necessario determinare, a partire dalle condizioni iniziali, il semiasse maggiore dell'ellissi. Questo può essere fatto a partire dall'energia meccanica totale. Sempre nel caso semplificato di orbite circolari si ha:
				H = (1/2) m (v^2) - 
            G (( m M)/r) = (1/2) m (ω^2) (r^2) - 
              G (( m M)/r) = (1/2) m ((G M)/(r^3)) (r^2) - 
                G (( m M)/r) = 
=-(1/2) 
                  G (( m M)/r)
Nel caso di orbite ellittiche questo risultato si generalizza a

					H = -(1/2) G (( m M)/a)
Da H possiamo calcolare a :

Pertanto per calcolare il periodo seguiamo il seguente schema:

		condizioni iniziali -> energia meccanica totale -> semiasse maggiore -> periodo
		
Notiamo che ciò ha senso se l'energia è negativa.

Perfezionare le istruzioni già viste in modo da costruire una funzione delle condizioni iniziali che restituisca l'orbita di un pianeta in unità MKS (G=6.67 (10^(-11)) N (m^2)/(Kg^2), M=2.01 (10^30) Kg) e che verifichi la conservazione dell'energia e del momento angolare.

Clear[Orbita]
\(Orbita[r0_, v0_] := 

      Module[{x, y, T, r, h, a, Fx, Fy, ax, ay, G = 6.67\ 10\^\(-11\), 
          M = 2.01\ 10\^30, En, 
          Ma}, 
h = \(1\/2\) 
              v0 . v0 - \(G\ M\)\/\@\(r0 . r0\); 
a = \(-\(1\/2\)\) \(G\ M\)\/h; 
T = \(\(2\ π\)\/\@\(G\ M\)\) \@a\^3; 
Print["Energia totale =", 
          h, "   Semiasse maggiore =", a, "   Periodo =", 
          T]; 
r[
            t_] = \@\(x[t]\^2 + y[t]\^2\); 
Fx = \(-\(\(G\ M\ m\)\/r[t]\^3\)\) 
            x[t]; 
Fy = \(-\(\(G\ M\ m\)\/r[t]\^3\)\) 
            y[t]; 
ax = Fx\/m; 
ay = 
          Fy\/m; 
{x, 
            y} = {x, 
              y} /. \(NDSolve[{\(x''\)[t] == ax, \(y''\)[t] == 
                    ay, x[0] == r0[\([1]\)], 
                  y[0] == 
                    r0[\([2]\)], 
\(x'\)[0] == 
                    v0[\([1]\)], \(y'\)[0] == v0[\([2]\)]}, {x, y}, {t, 
                  0, T}]\)[\([1]\)]; 

        ParametricPlot[{x[t], y[t]}, {t, 0, T}, 
          AspectRatio -> 1]; 
En[
            t_] := \(1\/2\) \((\(x'\)[t]\^2 + \(y'\)[t]\^2)\) - \(G\ M\)\/\@\(x[t]\^2 + y[t]\^2\); 
Plot[En[t], {t, 0, T}, 
          PlotRange -> {\(-10\^9\), 0}]; 
Ma[t_] := 
          x[t]\ \(y'\)[t]\  - \ y[t]\ \(x'\)[t]; 

        Plot[Ma[t], {t, 0, T}];
]\)

Esercizio. 
Si utilizzi la procedura ora sviluppata per					
       (a) Disegnare l'orbita della Terra. 
       (b) Determinare un set di condizioni iniziali per cui l'orbita non è chiusa

Clear[r0, v0]
r0 = {1.470994 (10^11), 0};
v0 = {0, 30286.5};
Orbita[r0, v0];

Clear[r0, v0]
r0 = {1.470994 (10^14), 0};
v0 = {0, 30286.5};
Orbita[r0, v0];