Laboratoriodifisicacomputazionale
			(Prova  d' esamedel18luglio2003)
Una particella di massa m si muove sotto l' 
              azione del potenziale
	\(U \((x)\) = \(\(U\_0\) \((\((x\/a)\)\^4 - \((x\/a)\)\^2)\) U\_0 = \(1  Joul/Kg  a = 1  m\)\)\)
	\(•\ \ \ \ Riportare\ in\ un\ grafico\ completo\ di\ etichette\ con\ unità\ di\ misure\ per\ gliassi\ il\ potenziale\ U \((x)\)\)
	\(•\ \ \ \ Scrivere\ in\ funzione\ delle\ condizioni\ iniziali\ le\ seguenti\ procedure\)
		(1) 
            calcolodel minimo e del massimo valore assuntoda x
		(2) calcolo del periodo del moto
		(3) calcolo dell' equazione oraria del moto
	\(•\ \ \ \ Riportare\ in\ un\ grafico\ completo\ di\ etichette\ con\ unità\ di\ misura\ per\ gli\ assi\)
	       \(l' equazione\ oraria\ del\ moto, 
          mettendo\ in\ evidenza\ le\ diverse\ situazioni\ \(\(fisiche\)\(.\)\)\)

Clear["Global`*"];
Off[General::"spell"];
Off[General::"spell1"];
Off[NIntegrate::"slwcon"]
U[x_] := (x^4) - (x^2);
a[x_] := -U'[x];
xlim[
          H_] = (Sqrt[z /. Solve[U[(Sqrt[z])] == H, z]]);
Bound[H_] := 
      Module[{x = xlim[H]}, 
        If[H < 0, 
          x, {-x[[2]], x[[2]]}]]
T[H_] := T[
          H] = 
Module[{f, a, b}, 
f[z_] := 
            If[H - U[z] == 0, 
              0. , (((2 (H - U[z])))^(-1/2))]; 
{a, b} = Bound[H]; 
Re[
            NIntegrate[
              f[z], {z, a, 
                b}]]
]
t[x_ /; NumberQ[x], 
        H_ /; NumberQ[H]] := t[x, 
          H] = 
Module[{a, b, 
            f}, 
{a, b} = 
            Bound[H]; 
f[z_] := 
            If[H - U[z] == 0, 
              0. , (((2 (H - U[z])))^(-1/2))]; 
Re[
            NIntegrate[f[z], {z, a, x}, 
              MaxRecursion -> 
                50]]
]
xs[H_ /; NumberQ[H]] := xs[
          H] = 
Module[{a, b, z1, z2, z, 
            n = 100, τ = T[H], ss = {}}, 
{a, b} = 
            Bound[H]; 
z1 = 
            a; 
Do[
z2 = 
              z1 + (b - a)/n; 
z0 = 
              Re[z /. FindRoot[
                    t[z, H] == ts, {z, z1, 
                      z2}]]; 
AppendTo[
              ss, {ts, z0}]; 
z1 = 
              z0, 
{ts, 
              0. , τ, τ/n}]; 
Interpolation[
            ss]
]
x[ts_, x0_, 
        v0_] := 
Module[{H = (1/2) (v0^2) + U[x0], 
            t0, τ, s, tt, y}, 
τ = 
            T[H]; 
y = xs[H]; 
s = 
            If[H < 0, Sign[x0], 1]; 
t0 = 
            ts + If[v0 > 0, t[s x0, H], 
                2  τ - t[s x0, H]]; 
tt = 
            Mod[t0, 2  τ]; 
s Which[
              0 < tt ≤ τ, y[tt], tt ≤ 2 τ, 
              y[2  τ - tt]]
] /;  
        NumberQ[ts + x0 + v0]
xsnum[
        H_ /; NumberQ[H]] := xsnum[
          H] = 
Module[{τ = T[H], Δt, 
            z1, v1 = 0. , z2, v2, ss, vv, 
            it}, 
z1 = Bound[
                H][[1]]; 
ss = {{0. , 
                z1}}; 
vv = {{0. , 
                v1}}; 
Δt = 
            0.0001  τ; 
Do[
{z2, 
                v2} = {(1/2) a[z1] Δ(t^2) + 
                  v1 Δt  + z1, 
                a[z1] Δt + v1}; 
AppendTo[
              ss, {it, z2}]; 
{z1, v1} = {z2, 
                v2}, {it, Δt, τ, Δt}
]; 
Interpolation[
            ss]
]
xnum[ts_, x0_, 
        v0_] := 
Module[{H = (1/2) (v0^2) + U[x0], 
            t0, τ, s, tt, y}, 
τ = 
            T[H]; 
y = xsnum[H]; 
s = 
            If[H < 0, Sign[x0], 1]; 
t0 = 
            ts + If[v0 > 0, t[s x0, H], 
                2  τ - t[s x0, H]]; 
tt = 
            Mod[t0, 2  τ]; 
s Which[
              0 < tt ≤ τ, y[tt], tt ≤ 2 τ, 
              y[2  τ - tt]]
] /;  
        NumberQ[ts + x0 + v0]

Plot[U[x], {x, -2, 2}, 
        AxesLabel -> {"x (metri)", "U (Joul/Kg)"}];

x0 = 0.9;
v0 = 0.5;
H = (1/2) (v0^2) + U[x0];
{(Subscript[x, 1]), (Subscript[x, 2])} = Bound[H];
τ = T[H];
\(\(Print[\*""\!\(x\_1\)="", x\_1, \*"" \!\(x\_2\)="", 
        x\_2, " T=", 2  τ];\)\)
Timing[xsnum[H]]
Timing[xs[H]]
Plot[{x[ts, x0, v0], xnum[ts, x0, v0]}, {ts, 0, 4  τ}, 
        AxesOrigin -> {0, 0}, PlotRange -> {0, 1}, 
        PlotStyle -> {RGBColor[0, 0, 1], RGBColor[0, 1, 0]}, 
        AxesLabel -> {"t (sec)", "x (metri)"}];

x0 = 0.9;
v0 = 1;
H = (1/2) (v0^2) + U[x0];
{(Subscript[x, 1]), (Subscript[x, 2])} = Bound[H];
τ = T[H];
\(\(Print[\*""\!\(x\_1\)="", x\_1, \*"" \!\(x\_2\)="", 
        x\_2, " T=", 2  τ];\)\)
Timing[xs[H]]
Timing[xsnum[H]]
Plot[{x[ts, x0, v0], xnum[ts, x0, v0]}, {ts, 0, 4  τ}, 
        AxesOrigin -> {0, 0}, PlotRange -> {(Subscript[x, 1]), (Subscript[x, 2])}, 
        PlotStyle -> {RGBColor[0, 0, 1], RGBColor[0, 1, 0]}, 
        AxesLabel -> {"t (sec)", "x (metri)"}];

x0 = -0.9;
v0 = 0.5;
H = (1/2) (v0^2) + U[x0];
{(Subscript[x, 1]), (Subscript[x, 2])} = Sort[Sign[x0] Bound[H]];
τ = T[H];
\(\(Print[\*""\!\(x\_1\)="", x\_1, \*"" \!\(x\_2\)="", 
        x\_2, " T=", 2  τ];\)\)
Timing[xs[H]]
Timing[xsnum[H]]
Plot[{x[ts, x0, v0], xnum[ts, x0, v0]}, {ts, 0, 4  τ}, 
        AxesOrigin -> {0, 0}, PlotRange -> {(Subscript[x, 1]), 0}, 
        PlotStyle -> {RGBColor[0, 0, 1], RGBColor[0, 1, 0]}, 
        AxesLabel -> {"t (sec)", "x (metri)"}];