Laboratorio di Fisica Computazionale

Capitolo 2. Calcolo Numerico

2.1 Operazioni aritmetiche

2.1 + 3.72

6.882/2

(54.8/4)

2^3
(2^3)

2.2  Costanti Matematiche

pi

N[Pi]

Pi
π

E
E

N[E]

2.3  Funzioni elementari

Sin[a]

Cos[h]

Round[1.8]

2.4 Calcoli standard

Sqrt[(1.1*10^-23)*(6.8*10^-2)/(1.4*10^-24)]

(Sqrt[((1.1*(10^(-23))*6.8*(10^(-2)))/(1.4*(10^(-24))))])

ScientificForm[0.0076]

2.5 Valori esatti ed approssimati

Precision[24]

345/678

(1/3)*3

N[(1/3)]*3

Precision[345/678]

x = 123/456  +  456/789  +  2 Sqrt[5]

Precision[x]

N[x]
x // N

Precision[N[x]]

Accuracy[N[x]]

Sin[3]
Sin[3. ]

x

Sqrt[x  +  4]
Sqrt[x  +  4. ]

x = (Sqrt[3.5])
Precision[x]

x = 2.5677

x = 
      111.34545645777777777777777777777777777777777774875858587584848

x

Precision[x]

Accuracy[x]

x = 1.0000001

x - 1

Cos[20 Degree]

Cos[20 Degree] // N

Cos[20] // N

2.6 Precisione arbitraria

a = (4 Pi/3) (2 m)^3

N[a]

N[a, 20]

N[(4 Pi/3) (2 m)^3,  200]

2.7 Matrici

Alcune funzioni che operano su matrici sono:
    c m    				prodotto di uno scalare c per una matrice m
    m.n   				prodotto della matrice m per la matrice n
    Inverse[m]
    MatrixPower[m,k]
    Det[m]
    Eigenvalues[m]
    Eigenvectors[m]
    MatrixForm[list]

Esempio 2.7.1

Determinare l'inversa della matrice
                                                       (({{16, 0, 0}, {0, 14, -6}, {0, -6, -2}}))

m  =  {{16,  0,  0}, {0,  
          14,  -6}, {0,  -6,  -2}}

MatrixForm[m]

minv  =  Inverse[m]

m . minv

MatrixForm[%]

Clear[m,  minv]

2.8 Doppie parentesi quadre

Esempio 2.8.1

Assegna il nome ourlist alla lista {1,3,5,7,9,11,13,15} estrai poi il secondo elemento e costruisci la lista formata dal primo, quarto e quinto elemento.

ourlist  =  {1,  uuu,  {2, 4},  7,  9,  11,  13,  
        15}

ourlist[[2]]

ourlist[[{1,  4,  5}]]

Esempio 2.8.2

Assegna il nome myMatrix alla matrice
								(({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}})) 
a) estrai la seconda riga;
b) estrai il secondo elemento della terza riga.

myMatrix  =  {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}

MatrixForm[myMatrix ]

myMatrix[[2]]

myMatrix[[1,  1]]

Clear[ourlist,  myMatrix]

2.9 Fit dei minimi quadrati

Esempio 2.9.1

Data la tabella (distanze-tempi)
			   			  t        d
						({{0, 0}, {1, 1.5}, {2, 6.0}, {3, 13.5}, {4, 24.0}, {5, 37.5}, {6, 54.0}, {7, 73.0}, {8, 96.0}, {9, 121.5}})
						
dove d è espresso in metri e t in secondi, determinare l'equazione che dà la distanza d in funzione del tempo t.

time  =  Table[i,  {i,  0,  9}]

\(distance\  = \ {0, \ 1.5, \ 6.0, \ 13.5, \ 24.0, \ 
        37.5, 
	\ 	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 54.0, \ 
        73.5, \ 96.0, \ 121.5}\)

val = {time, distance}

MatrixForm[val]

ListPlot[val]

data  =  Transpose[val]

MatrixForm[data]

ListPlot[data]

ListPlot[val]

fit = Fit[data,  {1,  t},  t]

fit = Fit[data,  {1,  t,  t^2,  t^3},  t]

Plot[fit, {t, 0, 9}]

Show[ListPlot[data], Plot[fit, {t, 0, 9}]]

2.10 Numeri Complessi

x = Sqrt[-7]  +  Log[2  +  8 I]

N[x]

Conjugate[N[x]]

Re[x]
Im[x]

Re[N[x]]
Im[N[x]]

2.11 Soluzione numerica di equazioni Polinomiali

NSolve[eqns,var] trova le soluzioni numeriche approssimate delle radici di una equazione polinomiale o di sistemi di equazioni polinomiali. In Mathematica una equazione si scrive nella forma lhs==rhs. Bisogna distinguere fra gli operatori "=" e "==":  "=" assegna un valore mentre "=="  è specifico per le equazioni.

x = 3

Clear[x]

x = 4

NSolve[x^5  +  4 x^3  + 3 x^2  +  2 x  - 10 ==  0,  
      x]

x

Clear[x]

NSolve[x^5  +  4 x^3  + 3 x^2  +  2 x  ==  10,  
      x]

x^5  +  4 x^3  + 3 x^2  +  2 x  ==  10

\(NSolve[{x^2\  - \ x\ y\  - \ y^2\  == \ 1, \ 
        x^3\  + \ 3\ x\ y\  - \ y^3\  == \ 
          9}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ {x, \ y}]\)

2.12 Integrazione Numerica

NIntegrate[f,{x,xmin,xmax}] da il valore approssimato dell'integrale

						(∫_xmin^xmax) f (x) dx
						
     Calcoliamo, per esempio, l'integrale
     						(∫_0^2) 
            sin[(Sqrt[1 + (x^4) + sin ((x^3))])] dx

NIntegrate[
      Sin[Sqrt[1  +  x^4  +  Sin[x^3]]],  {x,  0,  2}]

NIntegrate[f,{x,xmin,xmax},{y,ymin,ymax}] calcola il valore approssimato dell'integrale multidimensionale

						(∫_xmin^xmax)dx(∫_ymin^ymax) f (x, y) dy

Esempio 2.12.1

Un disco di raggio a e densità di carica σ  è posto in un piano xy con il suo centro nell'origine. Calcolare il potenziale elettrico nel punto P(1,2,5) (in unità di a).

\(\(\( (*\ 
      Generazione\ della\ figura\ *) \)\(
\)\(\(Unprotect[
        O];\)

    O = {0, 0, 0}; Px = {0, 1.5, 0}; Py = {1.5, 0, 0}; 
    Pz = {0, 0, 4};

    Q = {0.5, 0.6, 0}; P = {1, 2, 3}; Pt = P + {0, 0, 0.2}; 
    Qt = Q + {0.15, 0, 0}; r = \@\(Q . Q\)/2;

    \(fig1 = 
        ParametricPlot3D[{Cos[t], Sin[t], 0}, {t, 0, 2  π}, 
          DisplayFunction -> Identity];\)

    \(fig2 = 
        ParametricPlot3D[{r\ Cos[t], r\ Sin[t], 0}, {t, 0, 0.9}, 
          DisplayFunction -> Identity];\)

    \(g1 = 
        Show[fig1, 
          Graphics3D[{
Line[{O, 
                  Py}], 
Text[
                StyleForm["R", FontSize -> 6, 
                  FontWeight -> "Bold"], \((P + Q)\)/2 + {0.2, 0.2, 
                    0}], 
Text[
                StyleForm["r", FontSize -> 6, 
                  FontWeight -> "Bold"], {0.25, 0.6, 
                  0}], 
Text[
                StyleForm["φ", FontSize -> 6, 
                  FontWeight -> "Bold"], {0.25, 0.1, 
                  0}], 
Text[
                StyleForm["Q", FontSize -> 6, 
                  FontWeight -> "Bold"], 
                Qt], 
Text[
                StyleForm["P", FontSize -> 6, 
                  FontWeight -> "Bold"], 
                Pt], 
PointSize[0.02], Point[Q], 
              Point[P], 
Line[{O, Q}], 
              Line[{Q, P}]}], 
fig2, Axes -> None, 
          Boxed -> False, PlotRange -> All];\)

    \(g2 = 
        Show[Graphics[{
Text[
                StyleForm[\*""dV=\!\(1\/\(4  π\\ ϵ\_0\)\)(dq/R)"", FontSize -> 6, 
                  FontWeight -> "Bold"], {0, 
                  2}], 
Text[
                StyleForm["dq=σ dxdy=σ r drdφ", 
                  FontSize -> 6, FontWeight -> "Bold"], {0, 
                  1}], 
Text[
                StyleForm[\*""R=\!\(\@\(\((x - r\\ Cos \((φ)\))\)\^2 + \((y - r\\ Sin \((φ)\))\)\^2 + z\^2\)\)"", 
                  FontSize -> 6, FontWeight -> "Bold"], {0, 
                  0}]
}], Axes -> None, 
          DisplayFunction -> Identity, 
          PlotRange -> All];\)

    \(Figdisco = 
        Show[GraphicsArray[{g1, g2}, GraphicsSpacing -> 0.1], 
          DisplayFunction -> $DisplayFunction];\)

    \(Protect[O];\)\)\)\)

Il potenziale elettrico, dovuto al disco,  nel generico punto P(x,y,z) è:
				V(x,y,z)=(σ/(4  π (Subscript[ϵ, o])))(∫_0^a)\(dr \(∫\_0\%\(2  π\)r/√\((\((x - r\ Cos \((φ)\))\)\^2 + \((y - r\ Sin \((φ)\))\)\^2 + 
                      z\^2)\) dφ\)\)

((σ a)/(4  π(Subscript[ϵ, 0]))) 
      NIntegrate[
        r/Sqrt[(1  -  r Cos[phi])^2  +   (2  -  
                    r Sin[phi])^2  +  5^2], {r,  0,  1},  {phi,  
          0,  2  Pi}]

2.13 Soluzione Numerica di Equazioni Differenziali

Esempio 2.13.1

L'equazione del moto di un oscillatore armonico smorzato è
					(((d^2) x)/(dt^2)) + γ  (dx/dt) + (ω^2) x = 
        0

      
Posto ω=9 (s^(-2)), γ=0.5(s^(-1)), x(0)=1, x'(0)=0 determinare x(t)

sol = 
      NDSolve[{x''[t] + 0.5 x'[t] + 9 x[t] == 
            0, x'[0] == 0, x[0] == 1}, 
        x, {t, 0, 12}]

Plot[Evaluate[x[t] /. sol], {t, 0, 12}, 
      AxesLabel -> {"t (s)", "x (m)"}]

Clear[sol]

Esempio 2.13.2

L'intensità della forza gravitazionale dovuta al sole su un pianeta è
							F=((G M m)/(r^2))
					
Dove G, r, M, m sono la costante gravitazionele, la distanza del pianeta dal sole, la massa del sole e la massa del pianeta. La forza è diretta lungo la congiungente Sole-Pianeta.
Trattandosi di una forza centrale il moto avviene in un piano. Le equazioni del moto si possono quindi scrivere nella forma (l'origine del sistema concide con il centro del Sole):
\(\(\(									\)\(\(\(d\^2\) x\)\/dt\^2\)\)\)+((G M x)/((Sqrt[((x^2) + (y^2))])^3))=0
							(((d^2) y)/(dt^2))+((G M y)/((Sqrt[((x^2) + (y^2))])^3))=0

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)}];

Clear[orbit]

2.14 Iteratori

Gli iteratori sono dei comandi che dicono a Mathematica di ripetere delle istruzioni un certo numero di volte. Le notazioni possibili sono:
 {imax}							ripete imax volte
 {i,imax}						     ripete con i che assume i valori da 1 a imax
 {i,imin,imax}						     ripete con i che assume i valori da imin a imax
 {i,imin,imax,di}					     ripete con i che assume i valori da imin a imax incrementando i di di
  {i,imin,imax},  {j,jmin,jmax},... 			     ripete con i che assume i valori da imin a imax incrementando i di 1; per ogni i
                                                                                   j va da jmin a jmax con step di 1

Table[Sin[i] // N,  {i, 1,  13}]

Table[{i, (i^2)},  {i, 1,  13}]

(Subscript[x, k])

Clear[x]
Do[(Subscript[x, k]) = k, {k,  2,  4}]

(Subscript[x, 1])
(Subscript[x, 2])
(Subscript[x, 3])
(Subscript[x, 4])
(Subscript[x, 5])

2.15 Esercizi svolti

1. Calcolare l'espressione
\(\(\(							\)\(\([10\ \((\(10.8\ 10\^3\)\/300)\)\^\(1\/2\) + 
              4]\)\^\(1\/3\)\)\)\)

(((10 (((((10.8 (10^3))/300)))^((1/2))) + 4))^((1/3)))

2. Calcolare l'espressione
\(\(\(							\)\(\(\(10\^\(-24\)\ 10\^12\)\/10\^\(-14\)\) \@\(32000. \/2\^3\)\)\)\)

ScientificForm[(((10^(-24)) (10^12))/(10^(-14))) (Sqrt[(32000. /(2^3))]), 3]

3. Determinare la lunghezza del cateto a del triangolo

Clear[a];
a = ScientificForm[400.  Sin[35  °], 3]

4. Determinare la lunghezza del cateto b del triangolo

Clear[b];
b = ScientificForm[400.  Cos[35  °], 3]

5. Date le matrici 
				A=(({{1, 0, -1}, {2, 4, 7}, {5, 3, 0}}))                B=(({{6, 1}, {0, 4}, {-2, 3}}))
   determinare AB.

Clear[A, B];
A=(({{1, 0, -1}, {2, 4, 7}, {5, 3, 0}}));B=(({{6, 1}, {0, 4}, {-2, 3}}));
Print[
      MatrixForm[A . B]]

6. Data  la matrice C, determinare l'inversa

						C=(({{-2, 1, 3}, {0, -1, 1}, {1, 2, 0}}))

MatrixForm[Inverse[(({{-2, 1, 3}, {0, -1, 1}, {1, 2, 0}}))]]

7. Luapp Yllek, fisico marziano, ha scoperto un gas molto particolare. Egli ha misurato volume e pressione di questo gas ottenendo i seguenti dati (V,P), dove V è il volume misurato in litri marziani e p la pressione misurata in atmosfere marziane: (0.608,0.05), (0.430,0.10), (0.304,0.20), (0.248,0.30), (0.215,0.40), (0.192, 0.50). Determinare l'equazione della curva che meglio interpola il set di dati usando una combinazione di 1, (1/V) e (1/(V^2))

Clear[vp, f, V]; 
    vp = {{0.608, 0.05}, {0.430, 0.10}, {0.304, 0.20}, {0.248, 0.30}, {0.215, 
          0.40}, {0.192, 0.50}}; 
f[V_] = Fit[vp, {1, (1/V), (1/(V^2))}, V]

Clear[lp, fp]; 
    lp = ListPlot[vp, PlotStyle -> PointSize[0.02], 
        DisplayFunction -> Identity];
fp = Plot[f[x], {x, 0.05, 0.6}, DisplayFunction -> Identity]; 
    Show[lp, fp, DisplayFunction -> $DisplayFunction];

8. Usa il metodo dei minimi quadrati per trovare la migliore approssimazione lineare al seguente set di dati:
    (4,5), (6,8), (8,10), (9,12)

Clear[data];
data = {{4, 5}, {6, 8}, {8, 10}, {9, 12}};
f[x_] = Fit[data, {1, x}, x]

lp = 
        ListPlot[data, PlotStyle -> PointSize[0.02], 
          DisplayFunction -> Identity];
fp = Plot[f[x], {x, 0, 10}, DisplayFunction -> Identity]; 
    Show[lp, fp, DisplayFunction -> $DisplayFunction];

9. Determinare parte reale e parte immaginaria dell'espressione

						((4 + i)/(2 + 3 i))

Clear[a];
a = ((4 + I)/(2 + 3 I)); {Re[a], 
      Im[a]}

10. Esprimere il seguente numero complesso nella forma p + i q, dove p e q sono numeri reali con 15 cifre di accuratezza.
    			   \(\((3  π + 7  i)\)\ Cos[
            37  °] + \((2 + 
              8  i)\) e\^\(\(-3\) i + 2\)\)

Clear[a];
a = (3  π + 7  I) Cos[
            37  °] + (2 + 
              8  I) (E^(-3 I + 2)); SetAccuracy[a, 15]

11. Risolvere l'equazione

					(x^7) + (x^5) + 2 (x^3) + (x^2) + 1 = 0

NSolve[(x^7) + (x^5) + 2 (x^3) + (x^2) + 1 == 0, x]

12. Risolvere il sistema di equazioni
	
					(x^2) + x y  +  (y^2) = 1
					(x^3) + (x^2) y + x (y^2) + (y^3) = 4

NSolve[{(x^2) + x y  +  (y^2) == 1, 
        (x^3) + (x^2) y + x (y^2) + (y^3) == 4}, {x, y}]

13. Calcolare l'integrale
					(∫_0^1) ((x^3)/((e^x) - 
                1)) dx

NIntegrate[((x^3)/((E^x) - 1)), {x, 0, 1}]

14. Calcolare l'integrale
					\(TraditionalForm\`∫\_\(-∞\)\%∞\(\(\( 
                    H\_4\)(x)\)\^2\) \(e\^\(-x\^2\)\) dx\)
					
      dove (Subscript[H, n])(x) sono i polinomi di Hermite

Timing[
      NIntegrate[HermiteH[4, 
                x]^2 (E^(-(x^2))), {x, -∞, ∞}]]

Timing[
      Integrate[HermiteH[4, 
                x]^2 (E^(-(x^2))), {x, -∞, ∞}]]

15. Risolvere l'equazione differenziale
					
					(((d^2) y)/(d (x^2))) = 2x + y + 3((d y)/(d x))
					
       con y(2)=1 e y'(2)=-1. Calcola y(2.2) e riporta in un grafico y(x) per x compreso fra 2.0 e 2.3.

Clear[y]; 
    y = y /. NDSolve[{y''[x] == 2 x + y[x]  + 3 y'[x], 
              y[2] == 1, y'[2] == -1}, 
            y, {x, 0.5, 5}][[1]];
y[2.2]
Plot[y[x], {x, 0.5, 3}];

Clear[y]; 
    y = y /. DSolve[{y''[x] == 2 x + y[x]  + 3 y'[x], 
              y[2] == 1, y'[2] == -1}, y, 
            x][[1]];
y[x]
y[2.2]
Plot[y[x], {x, 0.5, 3}];

16. Risolvere l'equazione differenziale 

					(((d^2) y)/(d (x^2))) + (Sin^2)(x)((d y)/(d x)) + 3(y^2)= \(TraditionalForm\`\(e\^-\)\^\(x\^2\)\)
					
con y(0)=1 e y'(0)=0. Determinare y(1), y(2) e y(3) e disegnare il grafico di y(x) per x compreso fra 0 e 3.

Clear[x, y];
y = 
        y /. NDSolve[{y''[x] + Sin[x]^2 y'[x] + 
                    3  y[x]^2 == (E^(-(x^2))), 
                y[0] == 1, y'[0] == 0}, 
              y, {x, 0. , 3}][[1]];
{y[1. ], y[2. ], y[3. ]}
Plot[y[x], {x, 0. , 3}];

17. Risolvere l'equazione differenziale 

					(((d^2) y)/(d (x^2))) + 5 x((d y)/(d x)) -  (1-(x^3) ) y= 0
					
con y(0)=0 e y'(0)=1. Determinare y(1.8) e graficare y(x) per x compreso fra 0 e 4.

Clear[x, y];
y = 
        y /. NDSolve[{y''[x] + 
                    5  x y'[x] + (1 - (x^3)) y[x] == 0, 
                y[0] == 0, y'[0] == 1}, 
              y, {x, 0. , 4}][[1]];
y[1.8]
Plot[y[x], {x, 0. , 4}];

18. Si consideri il moto di un pianeta (G M=4(π^2)) dovuto alla forza di attrazione da parte del Sole

					F=(GmM/(r^2))
					
       (a) Disegnare la traiettoria per x(0)=1 y(0)=0, x'(0)=0 e y'(0)=8. L'orbita  è circolare? Spiegare
       (b) Disegnare la traiettoria per x(0)=1, y(0)=0, x'(0)=-2.5 e y'(0)=8. L'orbita è ellittica? Spiegare
       (c) Determinare un set di condizioni iniziali per cui l'orbita è una parabola
       (d) Determinare un set di condizioni iniziali per cui l'orbita è una iperbole

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

      Module[{GM = 4  π\^2, x, y, r, rr0 = Append[r0, 0], 
          vv0 = Append[v0, 0], e, l, ec, a, b, T, d, Fx, Fy, ax, ay}, 

        If[rr0 . vv0 == 0, 
          Print[\*""Poichè (Subscript[r, 0]) e (Subscript[v, 0]) sono ortogonali l'orbita può essere circolare.""]]; 

        l = \@\(Cross[rr0, \ vv0] . Cross[rr0, \ vv0]\); 

        e = \(1\/2\) vv0 . vv0 - GM\/\@\(rr0 . rr0\); 
d = l\^2\/GM; 

        ec = \@\(1 + \(2\ e\ l\^2\)\/\((GM)\)\^2\); 


        T = \(2  π \@\( rr0 . rr0\)\)\/\@\(vv0 . vv0\); 
        Print[T, ec, e, l]; 

        If[N[ec] == 0. , 
          Print[\*""essendo (Subscript[e, c])=0 l'orbita è circolare.""]]; 

        If[ec > 0 && \ ec < 1, a = \(-\(GM\/\(2  e\)\)\); 
          b = l \@\(-\(1\/\(2  e\)\)\); T = \(2  π\ a\ b\)\/l; 
          Print[\*""essendo (Subscript[e, c])="", 
            ec\ , " l'orbita è ellittica."]]; 

        If[N[ec] == 1. , 
          Print[\*""essendo (Subscript[e, c])=1 l'orbita è parabolica.""]]; 

        If[ec > 1, 
          Print[\*""essendo (Subscript[e, c])="", 
            ec\ , " l'orbita è iperbolica."]]; 

        r[t_] = \@\(x[t]\^2 + y[t]\^2\); 

        Fx = \(-\(\(GM\ m\)\/r[t]\^3\)\) x[t]; 

        Fy = \(-\(\(GM\ 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] == rr0[\([1]\)], 
                  y[0] == rr0[\([2]\)], \(x'\)[0] == 
                    vv0[\([1]\)], \(y'\)[0] == vv0[\([2]\)]}, {x, 
                  y}, {t, 0, T}]\)[\([1]\)]; 

        ParametricPlot[{x[t], y[t]}, {t, 0, T}, 
          AspectRatio -> 1]
]\)

Clear[r0, v0]
r0 = {1, 0};
v0 = {0, 2  π};
Orbita[r0, v0];

\(\(\(
\)\(Clear[r0, v0]

    \(r0 = {1, 0};\)

    \(v0 = {\(-2.5\), 2  π};\)

    \(Orbita[r0, v0];\)\)\)\)

Clear[r0, v0]
r0 = {1, 0};
v0 = {0, 2 (Sqrt[ 2]) π};
Orbita[r0, v0];

Clear[r0, v0]
r0 = {1, 0};
v0 = {0, 4  π};
Orbita[r0, v0];

19. L'equazione del moto dell'oscillatore di Van der Pol può essere scritta come

					(((d^2) x)/(d (t^2))) -ϵ (1 - (x^2))  ((d x)/(d t)) + x= 0
					
dove ϵ è un parametro positivo. Questa equazione del secondo ordine può essere ricondotta ad un sistema di due equazioni del primo ordine:
      					\(TraditionalForm\`\(\(
\)\(
\)\)\)((d x)/(d t)) = v
((d v)/(d t)) = ϵ (1 - (x^2)) 
            v  -  x
      Si consideri il caso con ϵ = 1. Risolvere numericamente il sistema di equazioni  e disegnare il diagramma di fase ( cioè v verso x) per 0 < t < 7 π   e per cinque sets di condizioni iniziali: v(0) = 0, e x(0) da 1 a 3 con un incremento di 0.5. Ripetere il calcolo per parecchie altre condizioni iniziali scelte a caso.
      Da questi diagrammi di fasi deriva il concetto di "ciclo limite".

\(\(\(
\)\(Clear[x, v, x0, t]

    \(v0 = 0;\)

    \(ϵ = 1;\)

    \(Do[{x[x0], 
            v[x0]} = {x, 
              v} /. \(NDSolve[{\(x'\)[t] == 
                    v[t], \(v'\)[
                      t] == ϵ \((1 - x[t]\^2)\) v[t] - x[t], 
                  x[0] == x0, v[0] == v0}, {x, v}, {t, 0. , 
                  7  π}]\)[\([1]\)]; 

        ParametricPlot[{\(x[x0]\)[t], \(v[x0]\)[t]}, {t, 0, 21}, 
          PlotRange -> {{\(-3\), 3}, {\(-3\), 3}}], {x0, 1, 3, 
          0.5}];\)\)\)\)

Clear[x, v, x0, t, i]
(Subscript[x, m]) = 5;
(Subscript[v, m]) = 2;
ϵ = 1;
\(\(Do[{x[i], 
            v[i]} = {x, 
              v} /. \(NDSolve[{\(x'\)[t] == 
                    v[t], \(v'\)[
                      t] == ϵ \((1 - x[t]\^2)\) v[t] - x[t], 
                  x[0] == \(x\_m\) Random[], 
                  v[0] == \(v\_m\) Random[]}, {x, v}, {t, 0. , 
                  7  π}]\)[\([1]\)]; 

        ParametricPlot[{\(x[i]\)[t], \(v[i]\)[t]}, {t, 0, 21}, 
          PlotRange -> {{\(-x\_m\), x\_m}, {\(-x\_m\), x\_m}}], {i, 1, 
          20}];\)\)

20. Un cerchio uniforme di raggio a e carica totale q è posto nel piano xy  con il centro nell'origine. Calcolare il potenziale elettrico nel punto P(2a,3a,5a) in unità di (q/(4  π(Subscript[ϵ, o]) a)) dove \(TraditionalForm\`ϵ\_o\ è\ la\ costante\ \(\(dielettrica\)\(.\)\)\)

Show[Figdisco]

Clear[r, θ]
A = {2 , 3, 5};
P = {r Cos[θ], r Sin[θ], 0};
(1/(4  π(Subscript[ϵ, 0]))) 
      (q/(a N[π])) NIntegrate[
        (1/(Sqrt[(P - A) . (P - A)])), {r, 0, 1}, {θ, 0, 
          2  π}]

2.16 Esercizi proposti

1. Fissare l'ipotenusa e un angolo adiacente  (in gradi) in un triangolo rettangolo e determinare i due cateti.

2. Costruire la tabella dei primi 30 numeri dispari.

3. Costruire una tabella di 10 integrali calcolati numericamente tra x=0 e x=1 di 10 funzioni f(x) diverse.

\(4. \ Costruire\ 20\ coppie\ di\ valori\ \ \((x\_i, 
        f \((x\_i)\))\)\ con\ f \((x)\)\ funzione\ a\ piacere\ e\ farne\ il\ grafico\ con\ ListPlot . \ 
        Fare\ anche\ il\ grafico\ di\ f \((x)\)\ e\ confrontarlo\ col\ \(\(precedente\)\(.\)\)\)

5.   Determinare numericamente le soluzioni di 2  (x^2) - 2  x + 
        8 = 0. 
\(Trovare\ le\ stesse\ soluzioni\ usando\ la\ formula\ esatta\ risolutiva\ delle\ equazioni\ di\ 2\^o\ grado\ \ trasformarlo\ in\ risultato\ numerico\ con\ \ 20\ cifre\ decimali\ \ e\ confrontarlo\ col\ \(\(precedente\)\(.\)\)\)

\(6. \ Fissare\ arbitrariamente\ i\ coefficienti\ e\ risolvere\ il\ sistema\ lineare\ nelle\ incognite\ x, 
    y : 
\(a\_\(1, 1\)\) x + \(a\_\(1, 2\)\) y = 
      b\_1\)
        (Subscript[a, 2, 1]) x + (Subscript[a, 2, 2]) y = (Subscript[b, 2])
\(usando\ NSolve . \ 
        Risolverlo\ \ anche\ col\ metodo\ di\ sostituzione\ e\ confrontare\ i\ \(\(risultati\)\(.\)\)\)