Laboratorio di Fisica Computazionale
Capitolo 7. Programmazione

7.1 Espressioni

Questa sezione discute i costituenti,  la manipolazione e la valutazione di espressioni.

7.1.1 Atomi

Molti oggetti di Mathematica sono composti da parti. Per esempio {a,b,c} è composto da quattro parti:

Table[Part[{a,b,c},i],{i,0,3}]

Ci sono anche oggetti, chiamati atomi, che non hanno parti :
		Tipo		Descrizione 			Esempio
		Symbol	Simbolo			x, E, Sin
		String		Stringa di caratteri		"Questa è una stringa"
		Integer		Intero				2
		Real		Numero reale approssimato	1.2345
				n.m dove n e m sono interi
		Razional	Numero razionale n/m dove	2/3
				n e m sono interi m≠0
		Complex	Numero Complesso nella	2.12 +3/5 i
				forma a + i b con a e b interi,
				reali o razionali.

Part[2/3,0]

Part[2.12+3/5 I,0]

La funzione AtomQ[object] ritorna True se object è un atomo False altrimenti.
La funzione Head[atom] dà il tipo di atom.

AtomQ[E]

Head[2.12+3/5I]

Part[(x^2)+1,0]

Head["Questa è una stringa"]

7.1.2 Rappresentazione interna

Mathematica è in grado di trattare diversi tipi di espressioni:
		Espressioni Algebriche				x+y
		Funzioni Matematiche				Sin[x]
		Equazioni 					y''[x]+y[x]+2x==0
		Operatori					Derivative[1][f]
		Assegnazioni					a=b
		Regole di trasformazioni			a->b
		Oggetti grafici					-Graphics-
		Liste						{a,b,c}

Mathematica rappresenta internamente le espressioni in modo uniforme: i costituenti sono atomi.  Vi è una regola per ottenere espressioni più complicate a partire da espressioni semplici : Se h, e1, e2,... sono espressioni anche h[e1,e2,... ] è una espressione, dove h è detta head e gli ei sono detti elementi dell'espressione. Una espressione può anche non contenere elementi, in tal caso essa ha la forma h[] .
In definitiva in Mathematica le espressioni sono gli atomi e tutti gli oggetti rapprentati internamente nella forma h[e1,e2,... ], detti espressioni normali. La funzione Head può anche essere usata per estrarre l' head  di una espressione.
La  Input Form di una espressione può differire dalla rappresentazione interna. Il comando FullForm[expr] mostra la rappresentazione interna di expr.

Head[x+y]

FullForm[x+y]

FullForm[y''[x]+y[x]+2x==0]

fig=Plot[x,{x,0,1}]

fig[[2]]

FullForm[fig]

Short[expr,n] dà una forma limitata a circa n righe di expr.

Short[FullForm[fig],4]

7.1.3 Manipolazione

Le posizioni delle parti in una espressione h[e1,e2,...] sono determinate in accordo alla forma interna: h è la parte 0 ei la parte i-sima. Se ei a sua volta è una espressione la posizione della j-sima parte è {i,j} e così via.
Le parti di una espressione che richiedono n indici per specificare la loro posizione sono dette di livello n. Il livello 0 si riferisce all'intera espressione.
Il comando Level[expr,levelspec] dà una lista di tutte le subespressioni di expr di livello specificato levelspec secondo il seguente schema
				n		i  livelli da 1 a n
				Infinity		tutti i livelli
				{n}		il livello n
				{n1,n2}	i livelli da n1 a n2
Se l'opzione Heads->True viene data come terzo argomento, Level ritorna anche le heads.

expr=1+x+2(x+3)^m

FullForm[expr]

Level[expr,{4}]

Level[expr,{4},Heads->True]

Table[Part[{a,b,c},i],{i,0,3}]

Tutte le funzioni che abbiamo introdotte per manipolare le liste possono essere usate per manipolare espressioni.
Vediamo qualche esempio

Select[1+x+2 (x^2)+3(x^3)+Sin[x],AtomQ]

FullForm[1+x+2 (x^2)+3(x^3)+Sin[x]]

expr=(1+x+2 (x^2)+3(x^3)+Sin[x])/.Plus->List

Select[expr,AtomQ]/.List->Plus

Select[1+x+2 (x^2)+3(x^3)+Sin[x],OddQ]

ReplacePart[x/y,1+x^2,{2,1}]

FullForm[x/y]

Position[x/y,y]

Insert[f[a,g[b],c,d],e,3]

Sort[f[d,c,b,a,a]]

Union[f[d,c,b,a,a]]

Partition[List[a,b,c,d,e],2]

Partition[f[a,b,c,d,e],2]

expr=f[a,h[b,h[c,d]],f[a,h[d]],h[e,f[g]],f[a,f[b,f[c,d]]]]

Flatten[expr]

expr1=expr/.f->List

Flatten[expr1]

expr

Flatten[expr,h]

expr=a+(b/c)+((d+(e/g))/(r+(s/t)))

FullForm[expr]

Position[expr,Power,{3}]

%/. 0->1

MapAt[f,expr,%]

MapAt[f,expr,Position[expr,Power,{3}]/.(0->1)]

7.1.4 Esercizi

1. Determina la rappresentazione interna delle seguenti espressioni
    (a) (Sqrt[(x^2)+8x+16])
    (b) (3/(x-1))+(2/(((x-1))^2))
    (c) Sin(x+y)
    (d) Sin(x+y)=Sin(x) Cos(y)+Cos(x) Sin(y)
    (e) (∫_0^(π/2))Sin xdx
    (f) ((a+i b)/(c-i d))

FullForm[(Sqrt[(x^2)+8x+16])]
FullForm[(3/(x-1))+(2/(((x-1))^2))]
FullForm[Sin[x+y]]
FullForm[Sin[x+y]==Sin[x] Cos[y]+Cos[x] Sin[y]]
FullForm[Hold[(∫_0^(π/2))Sin[ x]\[DifferentialD]x]]
FullForm[((a+I b)/(c-I d))]

2. Determina la rappresentazione interna delle seguenti espressioni
    (a) {a, {b, {c, d}}, {e, f}}
    (b) Integrate[(x^2),x]
    (c) Factor[(x^10) - (y^10)]
    (d) f[a]+g[c,d]/.f[x]->(x^3)
    (e) Sin /@ (a+b+c)
    (f) ((x^2)/(y^2))

FullForm[{a,{b,{c,d}},{e,f}}]
FullForm[Hold[Integrate[(x^2),x]]]
FullForm[Hold[Factor[(x^10) - (y^10)]]]
FullForm[f[a]+g[c,d]/.f[x]->(x^3)]
FullForm[Hold[Sin /@ (a+b+c)]]
FullForm[((x^2)/(y^2))]

3. Considera l'espressione
   				((x^2)/(x-1))==(1/(x-1))
   Ottieni una lista di sub-espressioni dal (1^o)al (3^o)livello

Level[((x^2)/(x-1))==(1/(x-1)),3]

4. Considera l'espressione
   			a+(b/c)+((d+e)/(1+(f/(1-g/h))))
   Estrai la sub-espressione g/h

expr=a+(b/c)+((d+e)/(1+(f/(1-g/h))));
-Apply[Part,Prepend[Position[expr,-(g/h)][[1]],expr]]

Position[expr,-(g/h)][[1]]

Prepend[Position[expr,-(g/h)][[1]],expr]

5. Nell'espressione
   			1+3x+3(x^2)+(x^3)+3 y+6 x y +3 (x^2)y +3(y^2)+3x (y^2)+(y^3)
   elimina tutti i termini contenenti y

6. Inserire b dopo ogni elemento nell'espressione

   			h[e1,e2,e3,e4,e5,e6,e7]

7. Dall' espressione
   			a+b/c+(d+e)/(f+g)
   ottenere l'espressione
   			a+b+1/c+d+e+1/(f+g)

8. Con l'istruzione Map elevare al quadrato entrambi i membri dell'espressione
   			Sqrt[1+x]==x^2+1

f[x_]:=x^2

exp=Sqrt[1+x]==x^2+1

Map[f,exp]

9. Applicare f ai numeratori dell'espressione
   			a+b/c+(d+e/f)/(r+s/t)

10. Trasformare in equazione la regola
				x^2+x->1

7.2 Pattern

La parola pattern si può tradurre con forma o modello. La potenza di Mathematica deriva dalla capacità di utilizzare il pattern-matching (concordanza delle forme). Abbiamo visto come si usa il pattern arg_ per definire le funzioni. Vedremo che il pattern _, chiamato  blank, è il pattern base e significa qualsiasi espressione. Si possono costruire altri pattern utilizzando blanks e nomi di pattern restrittivi. Possiamo così definire funzioni che accettano solo gli argomenti in accordo con un determinato pattern.

7.2.1 Blanks

Il pattern base è blank, che sta per qualsiasi espressione singola. La sua forma interna è "_" o Blank[]. Per esempio h[_] rappresenta la classe delle espressioni h[elem] dove elem può essere qualsiasi espressione singola

{h[],h[(x^2)],h[Sin[x]],h[3+4I],g[x],h[{b,c}],h[a,b,c,d]}/.h[_]->"MATCH"

Con il pattern _ possiamo costruire altri patterns. Per esempio il pattern {_,_} rappresenta la classe delle liste con due elementi

Position[{{},{a},{b,c},{c,{d,b}}},{_,_}]

Il pattern _^_ rappresenta la classe delle potenze.
La funzione Cases[{(Subscript[e, 1]),(Subscript[e, 2]),…}, … }, pattern, levelspec] dà la lista delle "(Subscript[e, i])" che soddisfano il pattern e che sono al livello levelspec
La funzione DeleteCases[{(Subscript[e, 1]),(Subscript[e, 2]),…}, pattern, levelspec] rimuove tutte le parti che soddisfano il pattern e che si trovano al livello  levelspec

Cases[{3,x,a^4,(Sqrt[1+x]),(a+b x)^(3n+1)},_^_]

FullForm[(Sqrt[1+x])]

DeleteCases[{3,x,a^4,(Sqrt[1+x]),(a+b x)^(3n+1)},_^_]

7.2.2 Pattern con un nome

Spesso è necessario assegnare un nome ad un pattern. I nomi dei pattern devono essere simboli; cioè Head[name] deve dare come risultato un simbolo.
Le forme bidimensionali non possono essere nomi di pattern perchè non sono simboli:

Head /@ {(x^y), (Subscript[x, y]), (Subscript[x, z]^y), (Underscript[x, y]), (SuperDagger[x]), (SuperPlus[x]), (SuperStar[x]), (OverHat[x]), (OverBar[x])}

Il pattern x_ sta per ogni singola espressione che viene chiamata x nel secondo membro di una definizione o regola di trasformazione. La sua forma interna è Pattern[x, Blank[]]

FullForm[x_]

(h[a] + h[b, c] + h[a, a]) h[d, e, f] /. h[x_, y_] -> x^y

(h[a] + h[b, c] + h[a, a]) h[d, e, f] /. h[x_, x_] -> x^x

myfunc[x_^n_] := {x, n}

Map[myfunc, {3, x, a^4, Sqrt[1 + x]}]

Clear[myfunc]

Per pattern più complessi x:pattern rappresenta qualsiasi espressione in accordo col pattern e che viene chiamata x al secondo membro di una definizione o regola di trasformazione

Sin[1 + a^2] /. h:Sin[x_ + y_] -> {h, x, y}

{Sin[a], Sin[a]} /. {h_, h:Sin[x_]} -> h

{Cos[a], Sin[a]} /. {h_, h:Sin[x_]} -> h

7.2.3 Patterns restrittivi

Un pattern rappresenta una classe di espressioni; vediamo come si possono costruire patterns che rappresentano sottoclassi di una classe di espresssioni.
Il pattern _h rappresenta qualsiasi singola espressione con head h. Il pattern x_h rappresenta qualsiasi singola espressione con head h chiamata x nel secondo membro di una definizione o regola di trasformazione.
Per esempio tramite le funzioni Union e  Cases, utilizzando pattern restrittivi, si possono estrarre le variabili di una espressione algebrica

myexpr=Expand[(x + y + z)^3]

Cases[myexpr, _Symbol, Infinity]//Union

Possiamo elevare al quadrato i numeri razionali di una lista

{x, "good", 3, 2/3, 5/7} /. x_Rational -> x^2

geomean[x_List] := Apply[Times, x]^(1/ Length[x])

geomean[{a1, a2, a3, a4}]

geomean[2, 3, 4]

Clear[geomean]

Consideriamo i polinomi di Legendre di ordine n 
					(Subscript[P, n])(x)=(1/((2^n)n!))(((d^n)((x^2)-1))/(dx^n))                        n>0 -1≤x≤1

legdPoly[n_Integer, x_] :=
	             	((1/(2^n n!)) D[(x^2 - 1)^n, {x, n}])//Simplify

Do[Print["n = ", i, "   ", legdPoly[i, x]], {i, 0, 3}]

legdPoly[2.5, x]

Si possono costruire pattern restrittivi utilizzando un test. Il pattern pattern?test rappresenta ogni espressione che soddisfa il pattern a cui l'applicazione di  test dà come risultato True.

RandomInteger[{1,100},10]

Count[%, _?EvenQ]

Clear[legdPoly]

legdPoly[n_Integer?NonNegative, x_] :=
	              ((1/(2^n n!)) D[(x^2 - 1)^n, {x, n}])//Simplify

legdPoly[2, x]

legdPoly[-2, x]

Infine il pattern pattern/;condition rappresenta qualsiasi espressione che soddisfa il pattern e per quale la valutazione di condition dia risultato True.

Cases[{a, {a, b}, {a, {b, c}}, 3, {3, 4}, {5, {6, 7}}},
            	x_List /; VectorQ[x, NumberQ], Infinity]

RandomInteger[{1,100},10]

Cases[%, x_Integer/;(x > 10 && x < 80)]

Clear[legdPoly]

legdPoly[n_Integer?NonNegative, x_ /; Abs[x] <= 1] :=
        ((1/(2^n n!)) D[(y^2 - 1)^n, {y, n}] /. y -> x)//Simplify

Plot[legdPoly[2,x],{x,-1,1},AxesLabel->{"x","legdPoly[2, x]"}]

Plot[legdPoly[2,x],{x,-2,2},AxesLabel->{"x","legdPoly[2, x]"}]

Abs[x] <= 1

legdPoly[1, x]

Come si vede ora legdPoly dà il risultato solo se Abs[x]≤1 . Si può risolvere questo problema definendo due funzioni che soddisfano pattern diversi

legdPoly[n_Integer?NonNegative,x_Symbol]:=((1/(2^n n!)) D[(x^2-1)^n,{x,n}])//Simplify

legdPoly[n_Integer?NonNegative, x_Real /; Abs[x] <= 1] :=
          ((1/(2^n n!)) D[(y^2 - 1)^n, {y, n}] /. y -> x)//Simplify

legdPoly[2, x]

legdPoly[2, 0.5]

newLegdP[n_Integer?NonNegative, x_Symbol] := LegendreP[n, x]

newLegdP[n_Integer?NonNegative, x_Real /; Abs[x] <= 1] :=
                                                                                                         LegendreP[n, x]

newLegdP[2, x]

newLegdP[2, 0.5]

Clear[legdPoly, newLegdP]

v[x_ /; Abs[x] > 1] := 0
v[x_ /; Abs[x] < 1] := -1

Plot[v[x],{x,-2.5`,2.5`},AxesLabel->{" x (a)"," V (Vo)"},PlotStyle->Thickness[0.01`]]

Clear[v]

v[x_ /; Abs[x] > a] := 0
v[x_ /; Abs[x] < a] := -1

v[0.5]

Clear[v]

{f[2], f[a], g[6], f[5]} /. f[x_] :> x /; NumberQ[x]

Notare che con :> l'espressione x /; NumberQ[x] viene valutata dopo che la regola di sostituzione è stata applicata.
Al contrario con -> la valutazione dell'espressione è immediata.

{f[2], f[a], g[6], f[5]} /. f[x_] -> x /; NumberQ[x]

func[x_, y_] := Sqrt[x - y] /; x > y

func[5, 4]

func[5, 5]

Clear[func]

func[x_, y_] = Sqrt[x - y] /; x > y

7.2.4 Equivalenze  Strutturali

Nel fare il confronto fra pattern ed espressione Mathematica guarda se vi è equivalenza strutturale tra la espressione ed il pattern. Espressioni matematiche equivalenti soddisfano patterns differenti se la loro struttura è diversa.
Per vedere se due espressioni hanno la stessa struttura, Mathematica guarda la forma interna delle espressioni

MatchQ[a^2 - b^2, x_^2 - y_^2]

MatchQ[Sin[a]^2 - Cos[b]^2, x_^2 - y_^2]

MatchQ[(a + b)(a - b), x_^2 - y_^2]

MatchQ[(a + b)(a - b), (x_ + y_)(x_ - y_)]

MatchQ[a - b, x_ - y_]

MatchQ[a - 2b, x_ - y_]

FullForm[a - 2b]

FullForm[x_ - y_]

7.2.5 Attributi

Alle funzioni si possono assegnare delle proprietà generali, dette attributi. Nella valutazione e nel pattern-matching Mathematica tiene conto degli attributi. La funzione Attributes[f] dà la lista degli attributi di f.  Ci sono due modi per assegnare un attributo ad una funzione
			Attributes[f]={attr1,attr2,...}       assegna gli attributi di f
			SetAttributes[f,attr]		       aggiunge attr agli attributi di f
L'espressione Attributes[f]={}  rimuove tutti gli attributi di f mentre ClearAttributes[f,attr] rimuove attr dagli attributi di f.
Gli attributi che influiscono sul pattern-matching sono OrderLess, Flat e OneIdentity. Qui esamineremo solo i primi due attributi.  Gli argomenti di una funzione con l'attributo OrderLess sono automaticamente ordinati secondo l'ordinamento standard. Per esempio se f ha tale attributo f[a,c,b] e f[b,a,c] e così via sono equivalenti a f[a,b,c]. In altre parole f è commutativa.
Addizioni e moltiplicazione sono commutative per cui le funzioni Plus e Times hanno l'attributo OrderLess. Vediamo come tale attributo influisce sul pattern-matching

{{a + b, a}, {a + b, b}} /. {x_ + y_, x_} :> {x, y}

Se una funzione f  ha l'attributo Flat allora le espressioni f[f[a],f[b,c]], f[f[a],b,f[f[c]]] sono tutte eqivalenti a f[a,b,c]. In altre parole f è associativa. Poichè addizioni e moltiplicazione sono associative le funzioni Plus e Times hanno anche l'attributo Flat. Vediamo con un esempio l'effetto combinato dei due attributi

(a + b + c) /.(x_ + y_) :> {x, y}

L'espressione soddisfa il pattern con x->a e y->b+c. Tuttavia, a causa degli attributi Flat e Orderless, ci sono molte altre possibilità: Mathematica ritorna il primo match che trova. Per vedere gli altri, alla regola possiamo unire una condizione che risulta sempre falsa, oppure utilizzare la funzione ReplaceList[expr,rules]:

(a + b + c) /.(x_ + y_) :> {x, y} /; Print[{{x}, {y}}]

ReplaceList[(a + b + c), (x_ + y_) -> {x, y}]

7.2.6 Defaults

Il pattern x_:v rappresenta una espressione che, se omessa, assume il valore v.
Il pattern x_h:v rappresenta una  espressione con head  h che, se omessa, assume il valore v

{a, b} /. {x_, y_:d} :> {x^2, y^2}

{a} /. {x_, y_:d} :> {x^2, y^2}

{a, 2} /. {x_, y_Integer:10} :> {x^2, y^2}

{a} /. {x_, y_Integer:10} :> {x^2, y^2}

{a, 2, 3} /. {x_, y_Integer:10} :> {x^2, y^2}

{a, b} /. {x_, y_Integer:10} :> {x^2, y^2}

Il pattern x_. rappresenta una espressione con valore default assegnato da Mathematica. Nel pattern x_+y_.  il valore default di y è 0 mentre nel pattern x_ y_. il valore default è 1.

a + b /. x_ + y_. :> x^2 + y^2

a /. x_ + y_. :> x^2 + y^2

a b /. x_ y_. :> x^2 + y^2

a /. x_ y_. :> x^2 + y^2

Esempio 7.2.1  Costruire un pattern che rappresenta qualsiasi differenza di due termini

MatchQ[a - 5, x_ - y_]

FullForm[x_ - y_]

FullForm[a - 5]

Cases[{a - Sqrt[b], a - 2b, a - b/3, a^2 - 5b^3, a - 5, 
	              Tan[x] - Sin[y], a + b, Tan[x] + Sin[x]},
	            x_ + n_?Negative y_.]

Esempio 7.2.2  Costruire un pattern che rappresenta qualsiasi espressione della forma (Subscript[a, n])(x^n) dove (Subscript[a, n]) è un numero, x una variabile e n un intero positivo.

MatchQ[3 x^2, a_ x_^n_]

MatchQ[x^2, a_ x_^n_]

MatchQ[2 x, a_ x_^n_]

Cases[{3 x^2, x^2, 2 x, x^1.5, 1/y, 3, a + b,
	              2 y^(-2), y^2 + 1, x/y^2, Sqrt[x], Tan[x]},
            	a_. x_Symbol^n_. /; NumberQ[a]&&Positive[n]&&IntegerQ[n]]

7.2.7 Patterns Ripetuti

Il pattern patt1|patt2|... è un pattern che ha una delle forme patti. Per esempio il pattern _Integer|_Symbol|_Rational rappresenta qualsiasi intero, variabile o numero razionale

Cases[{3, x, 2/3, 2 + 3 I, Tan[x], Sqrt[x], 1.5},
            	_Integer | _Symbol | _Rational]

Esempio 7.2.3: Costruire un pattern che rappresenti un monomio ( a (x^n)  n≥0)

Cases[{3 x^2, x^2, 2 x, x^1.5, 1/y, 3, a + b, 1.5,
	              2 y^(-2), y^2 + 1, x/y^2, Sqrt[x], Tan[x]},
     	      (a_. x_Symbol^n_. /; Positive[n]&&IntegerQ[n]&&NumberQ[a])
     	     |a_ /; NumberQ[a]]

p.. è un pattern che rappresenta una sequenza di una o più espressioni, ciascuna in accordo con p.

{a, a, a, a, b, b} /. {x_ .., y_ ..} :> {x, y}

7.2.8  Blanks Multipli

Il doppio blank rappresenta qualsiasi sequenza di una o più espression. L' input form è __ mentre la forma interna è BlankSequence[]. Il triplo blank ___ o BlankNullSequence[] rappresenta una sequenza di zero o più espressioni

myfunc[x_] := {x}

myfunc[a]

myfunc[a, b, c]

ourfunc[x__] := {x}

ourfunc[a]

ourfunc[a, b, c, d]

Il pattern __h rappresenta qualsiasi sequenza di una o più espressioni, ciascuna con head h. Quindi __List rappresenta una qualsiasi sequenza di una o più liste mentre {__List} rappresenta una lista di liste.

Cases[{{}, {a}, {a, b}, {{a}, {b, c}}}, {__List}]

{a, b, c, d} /.{x__, y__} :> "SMILE" /; Print[{{x}, {y}}]

Esempio 7.2.4: Calcolare la costante e con 60 cifre, determinare il numero di cifre tra il punto decimale ed una sequenza ripetuta di numeri e identificare la sequenza

N[E, 60]

ToString[%]

charactersOfE = Characters[%]

Print["digits  Sequence\n"];
Timing[
    charactersOfE /.
        {"2", ".", a___, x__/; Length[{x}] > 1, x__, ___} :>  "Ha!" /; Print[Length[{a}], "\t\t", x];
]

Clear[myfunc, ourfunc, charactersOfE]

7.2.9 Esercizi

1. Scrivi una funzione che calcola la media di una lista. Determina la media di una lista di n interi random nell'intervallo 0 1000, dove n stesso è un random fra 0 e 100.

mean[l_List]:=N[((l/.List->Plus)/(Length[l]))]
mean1[l_List]:=
Module[{i,k},k=Length[l];N[(1/k)(∑_(i=1)^k)l[[i]]]]

n=RandomInteger[{0,100}]
l=RandomInteger[{0,1000},n]
mean[l]
mean1[l]

FullForm[l]

2. Scrivi una funzione che espande solo l'argomento del logaritmo e quindi trasforma l'espressione 1-x+Log[(((a-b))^2)] in 1-x+Log[(a^2)-2 a b + (b^2)]

explog[expr_]:=expr/.Log[x__]:>Log[Expand[x]]
explog[1-x+Log[(((a-b))^2)]]

3. I polinomi associati di Legendre sono definiti dalla relazione
		 (Subscript[P, n]^m)(x)=(((-1))^m)(((1-(x^2)))^(m/2))(((d^m)(Subscript[P, n])(x))/(d (x^m)))  0≤m≤n
		 (Subscript[P, n]^m)(x)=(((-1))^m) (((n+m)!)/((n-m)!))(Subscript[P, n]^(-m))(x)                    -n≤m<0
Dove (Subscript[P, n])(x) è il polinomio d Legendre di ordine n e -1≤x≤1, n intero positivo.
Definisci una funzione che genera i polinomi associati di Legendre. Plot (Subscript[P, 2]^1)(x) e determina (Subscript[P, n]^m)(x) per n=0,1,2 e tutti i possibili valori di m. Verifica i risultati facendo uso della funzione intrinseca LegendreP.

Clear[AssLeg]
AssLeg[n_Integer/;n≥0,m_Integer/; m≥0,x_]:=Simplify[(((-1))^m) (((1-(x^2)))^(m/2)) Derivative[0,m][LegendreP][n,x]]/;m≤n
AssLeg[n_Integer/;n≥0,m_Integer/;m<0,x_]:=Simplify[(((((-1))^m) (m+n)! AssLeg[n,-m,x])/((n-m)!))]/;m≥-n

Plot[AssLeg[2,2,x],{x,-1,1}]

Table[AssLeg[n,m,x],{n,0,2},{m,-n,n}]
Simplify[Table[AssLeg[n,m,x]-LegendreP[n,m,x],{n,0,2},{m,-n,n}]]

Derivative[LegendreP[2,x]

7.3 Funzioni

7.3.1 Funzioni Pure

Abbiamo visto come in Mathematica si costruiscono funzioni che hanno un nome. Per evitare l'introduzione di molti simboli Mathematica permette la costruzione di funzioni anonime dette funzioni pure:
Function[var,body] è una funzione pura con un singolo parametro formale var. Function[var,body][arg] sostituisce arg al parametro formale var in body.
 Le due definizioni seguenti danno lo stesso risultato

f[y_]:=3 y

f[y]

Function[var,3 var][y]

Poichè il nome del parametro formale var è irrilevante, Mathematica permette l'uso di due forme ridotte: 
Function[body] e body&. In queste forme il parametro formale è  #.

Function[3#][y]

3#&[y]

Select[{1,a,x^2,3,5,1+x,7},#>4&]

test1[expr_]:=PolynomialQ[expr,x]
Select[1+x+2x^2+Sin[x],test1]

Select[1+x+2x^2+Sin[x],PolynomialQ[#,x]&]

q[x_]:=1/(1+x)

Nest[q,x,4]

Nest[1/(1+#)&,x,4]

Esempio 7.3.1  Dimostrare che
				▽(A.B)=A×(▽×B)+B×(▽×A)+(A.▽)B+(B.▽)A

Needs["VectorAnalysis`"]
SetCoordinates[Cartesian[x,y,z]];
A={AX[x,y,z],AY[x,y,z],AZ[x,y,z]};
B={BX[x,y,z],BY[x,y,z],BZ[x,y,z]};
Grad[A.B]==A\[Cross]Curl[B]+B\[Cross]Curl[A]+(A.{(∂_x)#1,(∂_y)#1,(∂_z)#1}&)/@B+(B.{(∂_x)#1,(∂_y)#1,(∂_z)#1}&)/@A

TableForm[Curl[B]]

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

func[1]

x

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

{r, s, x}

ClearAll[func, r, s, x]

7.3.2 Scelta della Definizione

In Mathematica è possibile dare ad un funzione più definizioni. Quando la funzione viene chiamata Mathematica parte dalla definizione più specifica e via via controlla quelle più generali. fino a quando la definizione è in accordo con la chiamata. Per esempio consideriamo le definizioni

ClearAll[f]
f[y_] := 3 y
f[x_Integer] := x^2
f[5] := 1/25

e la chiamata f[7]. Mathematica prima prova con f[5] che è la definizione più specifica poichè f[7] non concorda con f[5] Mathematica passa alla definizione f[x_Integer] che concorda con la chiamata. Quindi viene usata questa definizione

f[7]

?f

Spesso ci sono definizioni che hanno lo stesso ordine di generalità. In questo caso l' ordine è quello di immissione. Consideriamo per esempio le definizioni

ClearAll[f]
f[{x_, y_}] := x^3 - y^3
f[x_Real?Positive] := Sqrt[x]
f[x_Integer] := x^2

?f

La nozione di generalità di Mathematica è a volte incomprensibile. Si consideri per esempio le definizioni

ClearAll[f]
f[x_Integer] := x^2
f[x_Integer?Positive] := 1/x
f[x_Integer /; x > 50] := -x

ovviamente f[x_Integer /; x > 50]  è un caso particolare di f[x_Integer?Positive]  che a sua volta è un caso particolare di f[x_Integer]. Mathematica invece memorizza le definizioni secondo l'ordine di immissione.

?f

f[60]

ClearAll[f]
f[x_Integer?Positive] := 1/x
f[x_Integer /; x > 50] := -x
f[x_Integer] := x^2

?f

f[60]

ClearAll[f]

7.3.3 Funzioni Ricorsive e Programmazione Dinamica

Consideriamo, per esempio, la relazione di ricorrenza per i polinomi di Hermite
			(Subscript[H, n])(z) = 2(Subscript[H, n-1])(z) - 2(n-1)(Subscript[H, n-2])(z)		n>1
					(Subscript[H, 0])(z) = 1
					(Subscript[H, 1])(z) = 2 z
In questo caso la funzione e definita in termini di se stessa. Mathematica per l'uso di funzioni che chiamano se stesse.
Possiamo quindi calcolare i polinomi di Hermite con la seguente funzione ricorsiva:

hermite[n_Integer /; n > 1, z_] := Expand @ 
           (2 z hermite[n - 1, z] - 2 (n - 1) hermite[n - 2, z])

hermite[0, z_] := 1
hermite[1, z_] := 2 z

Timing[hermite[20, z];]

?hermite

Se l'ordine del polinomio n è grande il numero di chiamate alla funzione hermite[n,z] è molto alto. Per esempio  n=20 richiede 21891 chiamate.Se la funzione viene scritta in modo tale da ricordare in forma di regola o di definizione i valori calcolati, il numero di chiamate si riduce drasticamente. Il metodo di creare e memorizzare nuove regole o definizioni durante la valutazione è chiamato programmazione dinamica. Usando questa tecnica la definizione dei polinomi di Hermite diventa:

ClearAll[hermite]
hermite[n_Integer /; n > 1, z_] := hermite[n, z] = Expand @ 
                         (2 z hermite[n - 1, z] - 2 (n - 1) hermite[n - 2, z])
hermite[0, z_] := 1
hermite[1, z_] := 2 z

?hermite

hermite[8, z]

?hermite

ClearAll[hermite]

7.3.4 Iterazioni Funzionali

A volte si vuole applicare una funzione ripetutamente  ( un certo numero di volte o finchè un criterio è soddisfatto) ad una espressione. A questo scopo sono disponibile parecchie funzioni:

Nest[f,expr,n]				applica n volte f a expr e da il risulta finale
NestList[f,expr,n]			applica n volte f a expr e da la lista dei risultati {expr,f[expr],f[f[eexpr]],..}
FixedPoint[f,expr]			applica  f ripetutamente a expr finchè il risultato non cambia più
FixedPoint[f,expr,n]			Si ferma dopo n step. Il valore default di n è 65536
FixedPoint[f,expr,SameTest->comp]	Si ferma quando la funzione comp applicata a due successivi risultati da True 
FixedPoinList[f,expr]	
FixedPoinList[f,expr,n]
FixedPoinList[f,expr,SameTest->comp]

Per esempio il seguente comando genera nove termini di una progressione geometrica di ragione r

Plus @@ NestList[r #&, a, 8]

Possiamo implementare, utilizzando FixedPoint  l'algoritmo di Newton-Raphson per calcolare la radice quadrata di un numero r.  La formula da utilizzare è
				(Subscript[x, n+1])=(1/2)((Subscript[x, n])+(r/(Subscript[x, n])))
Partendo con (Subscript[x, 0])=3.0 calcoliamo (Sqrt[11])

r = 11;

FixedPoint[(# + r/#)/2&, 3.0]

% == Sqrt[11.0]

ClearAll[r]

Esempio 7.3.2  Consideriamo la mappa logistica data dall'equazione
					(Subscript[x, n+1])=f((Subscript[x, n]))
dove n=0,1,.... e 
					f(x)=μ x(1-x)
con 0≤x≤1 e 0≤μ≤4.
 (a) Mostrare che i punti x=0.6 e x=0 sono punti fissi stabili per μ=2.5 e μ=0.24 rispettivamente. 
 (b) Disegnare sullo stesso grafico (Subscript[x, n]) in funzione di n per due diversi punti iniziali (Subscript[x, 0])=0.500 e (Subscript[x, 0])=0.501

mu = 2.5;

f[x_] := mu x (1 - x)

(*0.6 == f[0.6]*)

Table[FixedPoint[f,RandomReal[]],{10}]

mu = 0.24;

0 == f[0]

FixedPoint[f, 0.612,SameTest->(Abs[#1-#2]<1*^-10&)]

FixedPointList[f, 0.612, 39]

%[[-1]]

SameQ[%[[-1]], %[[-2]]]

Table[Chop[FixedPoint[f,RandomReal[],SameTest->(Abs[#1-#2]<(1/(10^11))&)],(1/(10^10))],{10}]

Needs["ErrorBarPlots`"]

mu = 3.64;

ListPlot[{Transpose[{Range[0,44],NestList[f,0.5`,44]}],Transpose[{Range[0,44],NestList[f,0.501`,44]}]},Joined->True,PlotRange->{0,1},AxesLabel->{"n","xn"}]

Show[%,PlotRange->{{30,44},{0,1}}]

Nest[f, 0.500, 44]

Nest[f, 0.501, 44]

ClearAll[f, mu]

7.3.5 Protezione

Abbiamo visto che alle funzioni possiamo dare degli attributi. La maggior parte delle funzioni buit-in hanno l'attributo Protected: Se si cerca di aggiungere alle definizioni di una funzione con tale attributo un'altra definizione Mathematica non aggiunge la definizione e da un messaggio di errore. Per poter aggiungere definizioni bisogna prima rimuovere la protezione con il comando Unprotect[s1,s2,..] o ClearAttributes. Vediamo alcuni esempi.

Attributes[D]

D[u_ + v_, x_] := D[u, x] + D[v, x]

Clear[x]

D[x Sin[x], x]

Unprotect[D];

D[u_ v_, x_] := D[v, x] - v D[u, x]

Protect[D];

D[x Sin[x], x]

Unprotect[D];

D[u_ v_, x_] =.

Protect[D];

D[x Sin[x], x]

myfunc[x_] := x^2 + 1

Protect[myfunc];

Attributes[myfunc]

myfunc[x_, y_] := x^2 + y^2

Unprotect[myfunc];

myfunc[x_, y_] := x^2 + y^2

myfunc[a, b]

ClearAll[myfunc]

7.3.6 Upvalues and Downvalues

Abbiamo visto definizioni della forma
	f[args]=rhs
	f[args]:=rhs
Queste definizioni sono associate al simbolo f e sono dette downvalues per f.
Definizioni della forma
	f /: g[.....,f,....]=rhs 		o  (se g ha un solo argomento) g[f]^=rhs
	f /: g[.....,f,....]:=rhs		o  (se g ha un solo argomento) g[f]^:=rhs
		f /: g[.....,f[args],....]=rhs	o  (se g ha un solo argomento) g[f[args]]^=rhs
	f /: g[.....,f[args],....]:=rhs	o  (se g ha un solo argomento) g[f[args]]^:=rhs
sono anche associate ad f  ma sono dette upvalues per f. Nel caso in cui g ha più argomenti e si usa la seconda forma la definizione viene associata a tutti gli argomenti

Clear[f]
f[0] = 1;
f[x_, y_] := x^2 + y^3
f[g[x_], h[y_]] := p[x, y]

?f

f/: Re[f] = 0;
f/: Log[f[x_]] := q[x]
f/: g[f[x_], h[y_]] := w[x + y]

?f

ClearAll[a, b]

Abs[a + I b]

Im[a] ^= 0;
Positive[a] ^= True;

?a

Im[b] ^= 0;
Negative[b] ^= True;

Unprotect[Abs];

Abs[x_ /; (Im[x] == 0 && Positive[x])] := x
Abs[x_ /; (Im[x] == 0 && Negative[x])] := -x
Abs[u_^x_ /; Im[x] == 0] := Abs[u]^x
Abs[u_ v_] := Abs[u] Abs[v]
Abs[(x_ + y_. Complex[0, n_]) /;
  (MatchQ[Im /@ Level[{x, y}, {-1}], {(0)..}] &&
    FreeQ[{x, y},
                  w_^z_ /; (!IntegerQ[z] && !(Positive[w] === True))
               ]
  )
      ] := Sqrt[x^2 + (n y)^2]

FullForm[I]

Protect[Abs];

Abs[b]

Level[Sin[x], {-1}]

Flatten[{(3^(3/2))/(a^(2/3)) , (a^2/b^3) I}]

Im /@ Level[{(3^(3/2))/(a^(2/3)) , (a^2/b^3) I}, {-1}]

Abs[(3^(3/2))/(a^(2/3)) + (a^2/b^3) I]

Abs[((2/3) a - 3 b^2 a^(-7) I) b u^2/v^(2/3)]

Abs[(a^3 - (b^2/a^4) I)(u + a)^2/((2 + 3 a I)^5 v^3)]

7.3.7 Esercizi

1. Implementare il metodo di Newton-Raphson per trovare gli zeri di una funzione. Per determinare la radice di una funzione f(x), partire da un valore approssimato (Subscript[x, 0])e ottenere successivamente migliori approssimazioni con la formula 
				(Subscript[x, i+1])=(Subscript[x, i]) -((f((Subscript[x, i])))/(f'((Subscript[x, i]))))
Determinare i primi cinque zeri della funzione di Bessel del primo tipo (Subscript[J, 0])(x).

zero[f_,x0_]:=FixedPoint[#-((f[#])/(f'[#]))&,x0];
zeron[f_,x0_,n_]:=FixedPoint[#-((f[#])/(f'[#]))&,x0,n]

f[x_]:=(x^2)-1;
zero[f,0.9]
zero[((#)^2)-1&,0.8]

Plot[BesselJ[0,x],{x,0,15}]

Dynamic[MousePosition["Graphics"]]

f[x_]=BesselJ[0,x];
xstart={2.4,5.5,8.7,11.9,14.9}
Table[zero[f,xstart[[i]]],{i,Length[xstart]}]
Table[zeron[f,xstart[[i]],1],{i,Length[xstart]}]

2. Considerare la mappa a tenda definita dall'equazione
		(Subscript[x, n+1])=f((Subscript[x, n]))
dove n=0, 1, .... e
		f(x)=μ(1-2|x-(1/2)|)=2μ{({{x,, se 0≤x≤(1/2)}, {1-x, se (1/2)≤x≤1}})
con 0≤μ≤1.
a) Verificare che il punto x=0 è un punto fisso stabile quando il parametro μ<1/2.
b) Mostrare che vi sono due punti fissi quando μ >1/2 e nessuno di essi è stabile o attrattore.

7.4 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 π}]
x=y;

7.4.1 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 chiamata 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, "\n"];
                 p = 50; Print[p, "   ", q]]
p
q
Module[{q, p}, Print[q]; q = p^2; Print[q, "\n"];
                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=",v," u=",u]]
ClearAll[p, u]

7.4.2 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[cond1,val1,cond2,val2,....]	valuta in sequenza le condi e restituisce vali della prima che è vera.
Switch[expr,f1,v1,f2,v2,....]		valuta expr e la confronta con le fi, appena trova uguaglianza  restituisce vi

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"}]

Esempio: 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]

7.4.3 Iterazioni

Nel costruire iterazioni spesso si deve cambiare valore ad una variabile. Mathematica dispone delle seguenti notazioni per eseguire tali cambiamenti
 x++			aumenta x di 1 e restituisce il vecchio valore di x
++x			aumenta x di 1 e restituisce il nuovo valore di x
x--			diminuisce x di 1 e restituisce il vecchio valore di x
--x			diminuisce x di 1 e restituisce il nuovo valore di x
x+=dx			aggiunge dx a x e ritorna il nuovo valore
x-=dx			sottrae dx a x e ritorna il nuovo valore
x*=c			moltiplica  x per c e ritorna il nuovo valore
x/=c			divide  x per c e ritorna il nuovo valore
{x,y}={y,x}		scambia i valori di x e y
PrependTo[s,elem]	Inserisce s al primo posto in elem e aggiorna elem
AppendTo[s,elem]	Appende s ad elem e aggiorna elem

Le funzioni Do, While, For valutano procedure o espressioni ripetutamente. Il valore restituito da queste funzioni è Null.
La funzione Do valuta una espressione uno specificato numero di volte
Do[expr,{imax}]				valuta expr imax volte
Do[expr,{i,imax}]				valuta expr con la variabile i che assume successivamente i valori 1,2,..  imax
Do[expr,{i,imin,imax}]			il valore iniziale di i è imin
Do[expr,{i,imin,imax,di}]			usa incrementi di di
Do[expr,{i,imin,imax},{j,jmin,jmax}] 	a fissato i  fra imin e imax j assume i valori da jmin a jmax

x=1;
x++
x

x=1;
++x
x

sum=0;
Do[sum=sum+i;Print[sum],{i,1,5}]

ClearAll[f]
f[x_] := x^2
Do[f[i], {i, 3}]

Do[Print[f[i]], {i, 3}]

Do[Print["k = ", k, "    ",
	       Sum[1/(n^k), {n, 1, Infinity}]//N
	                ],
      {k, 2, 8, 2}
    ]

x = 0.500; Do[x = 3.64 x (1 - x), {99}]; x

Nest[(3.64 # (1 - #))&, 0.500, 99]

ClearAll[f, x]

Le funzioni While e For valutano una espressione o procedura ripetutamente finchè certe condizioni sono soddisfatte
While[condition,body]		valuta condition quindi body ripetutamente finchè condition diventa False
For[start,condition,incr,body]	esegue start quindi valuta ripetutamente body ed incr finchè condition diventa False

x = 1.0;
While[1 + 1/x =!= x, x = 1 + 1/x];
x

FixedPoint[(1 + 1/#)&, 1.0]

N[GoldenRatio]

ClearAll[x]

For[{n, v} = {1, {}}, Prime[n] < 40, n++,
         AppendTo[v, Prime[n]]];
v

ClearAll[n, v]

For[n = 1, n < 3, n++, Print[n]]

n

For[n = 1, n < 1, n++, Print[n]]

n

ClearAll[n]

7.4.4 Opzioni

Le funzioni possono avere due tipi di argomenti: argomenti posizionali e argomenti opzionali con un nome predefinito (opzioni). Per  esempio la funzione Plot ha numerose opzioni. Vediamo con un esempio come si definiscono le opzioni.
Come esempio scriviamo una procedura che calcola il potenziale ed il campo elettrico nel piano xy di un sistema di N cariche puntiformi disposte nello stesso piano.
			φ(x,y)=K (∑_(i=1)^N)((Subscript[q, i])/(Sqrt[(((x-(Subscript[x, i])))^2)+(((y-(Subscript[y, i])))^2)]))
			
			(Subscript[E, x])(x,y)=-((∂φ)/(∂x))
			
			(Subscript[E, y])(x,y)=-((∂φ)/(∂y))

Options[EField] ={SistemaDiMisura -> SI, angolo -> radianti};

EField[q:{{_,_,_}..},Optional[p:{_,_},{0,0}],opts___Rule]:=
Module[{φ,x,y,Ε,ℰ,θ,units,dir,n=Length[q]},
(* determina i valori delle opzioni  e li assegna a variabili locali*)
units=SistemaDiMisura/.{opts}/.Options[EField];dir=angolo/.{opts}/.Options[EField];		φ=(∑_(i=1)^n)((q[[i,1]])/(Sqrt[(((x-(q[[i,2]])))^2)+(((y-(q[[i,3]])))^2)]));
(Subscript[Ε, x])=N[-(∂_x)φ/.{x->(p[[1]]),y->(p[[2]])}];(Subscript[Ε, y])=N[-(∂_y)φ/.{x->(p[[1]]),y->(p[[2]])}];
ℰ=(Sqrt[(Subscript[Ε, x]^2)+(Subscript[Ε, y]^2)]);
θ=ArcTan[(Subscript[Ε, x]),(Subscript[Ε, y])];
φ=N[φ/.{x->(p[[1]]),y->(p[[2]])}];
Which[
units===SI,φ=φ 8.98755 (10^9) V;ℰ=ℰ 8.98755 (10^9) (V/m),
units===Gaussian,φ=φ statvolts;ℰ=ℰ(( statvolts)/cm),
True,Return["Il valore dell'opzione SistemaDiMisura deve essere SI o Gaussian."]
];
Print["Potenziale Elettrico = ",φ,"\n"];Print["Modulo del Campo Elettrico  = ",ℰ,"\n\n"];Which[
dir===radianti,θ=N[θ] rad,
dir===gradi,θ=N[(θ/Degree)] Gradi,
True,Return["Il valore dell'opzione angolo deve essere radianti o gradi."]
];
Print["Angolo che il campo elettrico forma con l'asse x  = ",θ];
]

EField[{{500, -10, 30}, {-200, 30, 20}},  SistemaDiMisura -> Gaussian, angolo -> gradi]

EField[Table[{2 10^-6, 0, i}, {i, -50, 50}], {100, 0}]

SetOptions[EField,SistemaDiMisura -> Gaussian];

EField[{{500, -10, 30}, {-200, 30, 20}},angolo -> gradi]

ClearAll[electrostaticField]

7.5 Problemi di Meccanica

7.5.1 Costruzione di una procedura per l'integrazione numerica delle equazioni del moto

Consideriamo il caso del moto in una dimensione di un corpo soggetto ad una forza elastica:

                (( (d^2)x)/(d (t^2)))=(F/m)    ;       F=-k x    
Costruiamo una sequenza equivalente a NDSolve.

Dati e condizioni iniziali:

Clear[f,m,t0,t1,x0,v0]
f[x_,v_,t_]=-2*x;
m=.5;
t0=0.;
t1=10.;
x0=0.;
v0=.5;

Troviamo la soluzione con NDSolve:

Clear[x,t]
sol = NDSolve[{x''[t] - ((f[x[t],x'[t],t])/m) ==0, x'[0]==v0, x[0]==x0}, x, {t, t0, t1}]

Flatten[sol]

Clear[xes];
xes[t_]:=x[t]/.Flatten[sol]

plotxes=Plot[xes[t],{t,t0,t1},PlotRange->All,PlotStyle->RGBColor[1,0,0]]

Il programma equivalente si ottiene iterando la sequenza: 

(( Δ x(t))/(Δ t))≃(( d x(t))/(d t))=v(t)    ;(( Δ v(t))/(Δ t))≃(( d v(t))/(d t))=((F(x,v,t))/m)
cioè
x(t+dt)= x(t) +v(t) dt   ;  v(t+dt)= v(t) +((F(x,v,t))/m) dt   

con dt fissato e sufficientemente piccolo:

npoint=100;
dt=((t1-t0)/npoint);
Clear[xx,vv];
xx[0]=x0;
vv[0]=v0;

Do[
xx[k+1]=xx[k]+vv[k]*dt;
vv[k+1]=vv[k]+((f[xx[k],vv[k],t])/m)*dt,
{k,0,npoint}]

xvec=Table[{t0+k*dt,xx[k]},{k,0,npoint}];

plotxvec=ListPlot[xvec,PlotRange->All]

Show[plotxes,plotxvec,PlotRange->All]

Al crescere di npoint, cioè al diminuire di dt, la soluzione approssima sempre meglio quella esatta. 
Abbiamo ottenuto una soluzione xvec per punti. Interpoliamola:

Clear[xap];
xap[100]=Interpolation[xvec];

xvec[3.67]

xap[100][3.67]

Plot[xap[100][t],{t,t0,t1},PlotRange->All]

npoint va scelto come compromesso tra velocità e precisione. Perfezioniamo il programma: partiamo da un npoint piccolo e lo raddoppiamo ogni volta finché l' "errore medio" definito dalla media del  valore assoluto della differenza di due soluzioni successive

(((∫_to^t1) | (Subscript[xap, 2 npoint])(t)-(Subscript[xap, npoint])(t)|d t)/(t1-t0))
non è uguale almeno a un valore prefissato ϵ .

npoint=200;
dt=((t1-t0)/npoint);
Clear[xx,vv];
xx[0]=x0;
vv[0]=v0;

Do[
xx[k+1]=xx[k]+vv[k]*dt;
vv[k+1]=vv[k]+((f[xx[k],vv[k],t])/m)*dt,
{k,0,npoint}]

xvec=Table[{t0+k*dt,xx[k]},{k,0,npoint}];

xap[200]=Interpolation[xvec]

Plot[{xap[100][t],xap[200][t]},{t,t0,t1},PlotRange->All]

Plot[Abs[xap[100][t]-xap[200][t]],{t,t0,t1},PlotRange->All]

Inseriamo per questo il Do precedente in un nuovo Do nel quale si calcola l'errore e che viene interrotto dal controllo If appena l'errore prefissato viene raggiunto:

ϵ=.01;(* errore richiesto *)
nini=10; (* npoint iniziale *)
npoint=(nini/2);
Clear[xapold];
xapold[t_]=0;
Do[
npoint=2*npoint;
dt=((t1-t0)/npoint);
Clear[xx,vv];
xx[0]=x0;
vv[0]=v0;

Do[
xx[k+1]=xx[k]+vv[k]*dt;
vv[k+1]=vv[k]+((f[xx[k],vv[k],t])/m)*dt,
{k,0,npoint}];

xvec=Table[{t0+k*dt,xx[k]},{k,0,npoint}];
Clear[xapnew];
xapnew=Interpolation[xvec];
err=(1/(t1-t0))*NIntegrate[Abs[xapnew[tt]-xapold[tt]],{tt,t0,t1},MaxRecursion->200];
Print[npoint,"    ",err];
Clear[xapold];
xapold[t_]=xapnew[t];
If[err≤ϵ,Break[]],
{iter,1,20}];

Clear[xap];
xap[t_]=xapold[t];

plotxap=Plot[xap[t],{t,t0,t1},PlotRange->All]

Show[plotxes,plotxap]

Esercizio: Inserire le istruzioni in una procedura che calcola la legge oraria. Usare il costrutto:
x[k_, m_, x0_, v0_, ϵ_]:=Module[{ ... }, ....]

Esercizio: aumentare gradualmente il valore della costante elastica: valutare come cambia l'errore e spiegare.

Esercizio: ripetere l'esempio aggiungendo una forza di attrito.

7.5.2 Problema: Moto di un pianeta

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. 


(afelio e perielio misurati in milioni di Chilometri)
In forma vettoriale la II Legge della Dinamica si scrive: 
			(OverBar[F])=-G((m M)/(r^3))(OverBar[r])=m(OverBar[a])=m(((d^2)(OverBar[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

Ovviamente per disegnare una e una sola orbita occorre conoscere il periodo. Lo ricaviamo nel caso più semplice di orbite circolari. In questo caso l'accelerazione è tutta accelerazione 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/m)=(1/2) (v^2)-G(M/r)=
=(1/2)(ω^2)(r^2)-G(M/r)=(1/2) ((G M)/(r^3)) (r^2)-G(( M)/r)=
=-(1/2)G(( M)/r)
Nel caso di orbite ellittiche questo risultato si generalizza a

					(H/m)=-(1/2)G(( M)/a)
Da H possiamo calcolare a :
			a=-(1/2)G(M/(H/m))
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[{r,h,a,Fx,Fy,ax,ay,G=(6.67/(10^11)),M=2.01 (10^30),En,Ma},h=((v0.v0)/2)-((G M)/(Sqrt[r0.r0]));a=-((G M)/(2 h));If[a>0,T=(((2 π) (Sqrt[(a^3)]))/(Sqrt[G M]));Print["Energia totale =",h,"   Semiasse maggiore =",a,"   Periodo =",T],T=((2 (2 π) (Sqrt[-(a^3)]))/(Sqrt[G M]));Print["Energia totale =",h,"Orbita aperta"]];r[t_]=(Sqrt[((x[t])^2)+((y[t])^2)]);Fx=-(((G M m) x[t])/((r[t])^3));Fy=-(((G M m) y[t])/((r[t])^3));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]];
Print[ParametricPlot[{x[t],y[t]},{t,0,T},PlotRange->All]];
En[t_]:=(1/2) (((x'[t])^2)+((y'[t])^2))-((G M)/(Sqrt[((x[t])^2)+((y[t])^2)]));Print[Plot[En[t],{t,0,T},PlotRange->Sort[{2 En[0],0}]]];
Ma[t_]:=x[t] y'[t]-y[t] x'[t];
Print[Plot[Ma[t],{t,0,T},PlotRange->Sort[{0,2 Ma[0]}]]];{x,y}]

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]
Clear[x,y,T]
r0={1.470994 (10^11),0};
v0={0,30286.5};
{x,y}=Orbita[r0,v0];

Clear[r0,v0]
Clear[x,y,T]
r0={1.470994 (10^12),0};
v0={0,30286.5};
Orbita[r0,v0]

7.5.3 Problema: Moto di due corpi soggetti a forze generiche.

Le equazioni del moto in 2 dimensioni di 2 corpi soggetti ad una forza reciproca (Subscript[F, 12])
piu' una forza esterna su ciascuno , (Subscript[F, 1]) ed (Subscript[F, 2]) , sono:

                (( (d^2)(Subscript[x, 1]))/(d (t^2)))=(((Subscript[F, 1]) + (Subscript[F, 12]))/(Subscript[m, 1]))    ;      (( (d^2)(Subscript[x, 2]))/(d (t^2)))=(((Subscript[F, 2]) - (Subscript[F, 12]))/(Subscript[m, 2]))

Risolviamole    ed  animiamo la soluzione nel caso di due corpi soggetti ad una attrazione gravitazionale reciproca che si muovono nel campo gravitazionale di un terzo corpo ( p. es. un sistema pianeta-satellite nel campo di una stella fissa ).

Clear[f,TwoBody,Filmato,M,m,m1,m2,r,r10,r20,v10,v20,t]
f[M_,m_,r:{_,_},Optional[v:{_,_},{0,0}],t___]=-((m M)/(((r.r))^(3/2)))r;
TwoBody[m1_,m2_,M_,r10:{_,_},r20:{_,_},v10:{_,_},v20:{_,_},t0_,t1_]:=
TwoBody[m1,m2,M,r10,r20,v10,v20,t0,t1]=
Module[{eq,r1,r2,x1,y1,x2,y2},
r1[t_]={x1[t],y1[t]};
r2[t_]={x2[t],y2[t]};
eq=Flatten[Map[{#[[1,1]]==#[[2,1]],#[[1,2]]==#[[2,2]]}&,
{m1 r1''[t]==f[M,m1,r1[t]]+f[m1,m2,r1[t]-r2[t]],
m2 r2''[t]==f[M,m2,r2[t]]+f[m1,m2,r2[t]-r1[t]],
r1[0]==r10,r2[0]==r20,r1'[0]==v10,r2'[0]==v20}]];
{x1,y1,x2,y2}={x1,y1,x2,y2}/.NDSolve[eq, {x1,y1,x2,y2}, {t, t0, t1},MaxStepSize->.1,StartingStepSize->.1,MaxSteps->100000][[1]]
];
Filmato[t0_,t1_,np_,d_]:=
Module[{dt=((t1-t0)/np)},
Table[
Show[Graphics[{RGBColor[0,0,0],PointSize[.09],Point[{0,0}]}],Graphics[{RGBColor[0,0,0],PointSize[.06],Point[{x1[t],y1[t]}]}],Graphics[{RGBColor[1,0,0],PointSize[.03],Point[{x2[t],y2[t]}]}],Axes->True,PlotRange->{{-d,d},{-d,d}},AspectRatio->Automatic],
{t,t0,t1,dt}]
];

m1=1;
m2=0.1`;
M=5;
r10={0,30};
r20={0,39};
v10={0.38`,0};
v20={0.05`,0};
t0=0.`;
t1=3000.`;
{x1,y1,x2,y2}=TwoBody[m1,m2,M,r10,r20,v10,v20,t0,t1];
Plot[{x1[t],x2[t]},{t,t0,t1},PlotRange->All]
Plot[{y1[t],y2[t]},{t,t0,t1},PlotRange->All]

Per animare la soluzione generiamo una sequenza di istantanee:

t0=0.;
t1=3000.;
np=500;
d=50;
fig=Filmato[t0,t1,np,d]
ListAnimate[fig]

Export["N:\\MyDir\\Didattica\\Fisica Computazionale\\Mathematica 6.0\\duecorpi.gif",fig,"GIF"]

7.5.4 Problema: Urto anelastico.

Il programma è analogo al precedente. Abbiamo ora solo due corpi sferici di raggio (Subscript[r, 1]) ed (Subscript[r, 2]): il primo si conficca nel secondo perdendo energia a causa di una forza di attrito per esempio proporzionale alla velocità relativa del primo corpo rispetto al secondo:

(Overscript[F, ->]) = - k ( (Overscript[(Subscript[v, 1]), ->]) - (Overscript[(Subscript[v, 2]), ->]) )      con    k > 0       se | (Overscript[(Subscript[x, 1]), ->])-(Overscript[(Subscript[x, 2]), ->]) | < (Subscript[r, 1])+ (Subscript[r, 2])

(Overscript[F, ->]) = 0      altrimenti.

Tutte le altre forze sono nulle.

Clear[f,Fa,Urto,FilmUrto,M,m,m1,m2,r,r10,r20,v10,v20,t]
Urto[f_,m1_,m2_,r0_,r10:{_,_},r20:{_,_},v10:{_,_},v20:{_,_},t0_,t1_]:=
Urto[m1,m2,r0,r10,r20,v10,v20,t0,t1]=
Module[{eq,r1,r2,x1,y1,x2,y2},
r1[t_]={x1[t],y1[t]};
r2[t_]={x2[t],y2[t]};
eq=Flatten[Map[{#[[1,1]]==#[[2,1]],#[[1,2]]==#[[2,2]]}&,{m1 r1''[t]==f[r1[t]-r2[t],r1'[t]-r2'[t],r0],
m2 r2''[t]==-f[r1[t]-r2[t],r1'[t]-r2'[t],r0],
r1[0]==r10,r2[0]==r20,r1'[0]==v10,r2'[0]==v20}]];
{x1,y1,x2,y2}={x1,y1,x2,y2}/.NDSolve[eq, {x1,y1,x2,y2}, {t, t0, t1},MaxStepSize->.1,StartingStepSize->.1,MaxSteps->100000][[1]]
];
FilmUrto[t0_,t1_,np_,d_]:=
Module[{dt=((t1-t0)/np)},
Animate[
Show[
Graphics[{RGBColor[1,0,0],PointSize[.06],Point[{x2[t],y2[t]}]}],Graphics[{RGBColor[0,0,0],PointSize[.02],Point[{x1[t],y1[t]}]}],Axes->True,PlotRange->d,AspectRatio->Automatic],
{t,t0,t1,dt}]
];

Clear[x1,x2,y1,y2,Fa]
Fa[r:{_,_},v:{_,_},r0_]={If[(Sqrt[r.r])≤r0,-0.2v[[1]],0],If[(Sqrt[r.r])≤r0,- 0.2v[[2]],0]};
m2=1.`;
m1=7.`;
r1=1.`;
r2=5.`;
t0=0.`;
t1=150.`;
r10={0,0};
r20={50.`,0};
v10={1.`,0};
v20={0,0};
{x1,y1,x2,y2}=Urto[Fa,m1,m2,r1+r2,r10,r20,v10,v20,t0,t1]
Plot[{x1[t],x2[t]},{t,t0,t1},PlotRange->All]
Plot[{y1[t],y2[t]},{t,t0,t1},PlotRange->All]

np=100;
d={{-20,130},{-20,20}};
FilmUrto[t0,t1,np,d]

Esercizio: Mostrare graficamente che la quantità di moto del sistema si conserva.
Mostrare graficamente che l'energia cinetica del sistema non si conserva.

7.5.5 Problema: Urto elastico.

Abbiamo due corpi di raggio r1 ed r2 : la forza tra essi e' repulsiva di modulo molto grande e diretta lungo la congiungente i loro centri (Overscript[(Subscript[x, 1]), ->]) e (Overscript[(Subscript[x, 2]), ->]) non appena si toccano ; nulla altrimenti. Detta r la distanza tra i centri possiamo rappresentarla

(Overscript[F, ->]) = - k ( (Overscript[(Subscript[x, 2]), ->]) - (Overscript[(Subscript[x, 1]), ->]) )      con    k > 0       se | (Overscript[(Subscript[x, 1]), ->])-(Overscript[(Subscript[x, 2]), ->]) | < (Subscript[r, 1])+ (Subscript[r, 2])

(Overscript[F, ->]) = 0      altrimenti.

Clear[x1,x2,y1,y2,Fe]
Fe[r:{_,_},v:{_,_},r0_]=Module[{K=(10^4)},{If[(Sqrt[r.r])≤r0,-K r[[1]],0],If[(Sqrt[r.r])≤r0,-K r[[2]],0]}];
m1=1.`;
m2=1.`;
r1=2.`;
r2=2.`;
t0=0.`;
t1=2.`;
r10={6,9};
r20={18,7};
v10={20.3`,0.9`};
v20={0.`,0.`};
{x1,y1,x2,y2}=Urto[Fe,m1,m2,r1+r2,r10,r20,v10,v20,t0,t1]
Plot[{x1[t],x2[t]},{t,t0,t1},PlotRange->All]
Plot[{y1[t],y2[t]},{t,t0,t1},PlotRange->All]

np=100;
d={{-20,40},{-20,20}};
FilmUrto[t0,t1,np,d]

Esercizio: dimostrare graficamente che la quantita' di moto e l'energia cinetica del sistema si conservano.

7.5.6 Problema: Moto di un razzo.

Un razzo di massa (Subscript[M, r]) contenente una massa (Subscript[M, c]) di carburante  viene lanciato da una rampa lunga d,inclinata di un angolo θ rispetto all'orizzontale. Il combustibile viene espulso con una velocità relativa al razzo v e la rapidità di esplusione del combustibile è k. Trascurando la resistenza dell'aria e supponendo costante la forza gravitazionale:
    - calcolare la gittata
	- riportare la traiettoria  in un grafico completo di etichette con unità di misura per gli assi 
	-  animare il moto del razzo.
Utilizzare i seguenti dati: (Subscript[M, r])=700 Kg,(Subscript[M, c]) =2100 kg,k=90 kg/s,d=10m ,θ=(π/4),v=2.4 (10^3) m/s .
Suggerimento: considerare le tre fasi del moto:
1) moto sulla rampa,
2)moto soggetto alla forza peso e alla spinta del motore,
3)moto,una volta esaurito il carburante,soggetto unicamente alla forza peso.

Clear[time,Moto,Filmato]
time[t_]:=
Module[{i,ore,minuti,secondi,decimi,fm},
	fm[i_]:=If[StringLength[ToString[i]]==2,ToString[i],"0"<>ToString[i]];
	ore=IntegerPart[t/3600];
	minuti=IntegerPart[(t-3600 ore)/60];
	secondi=IntegerPart[t-3600 ore-60 minuti];
	decimi=IntegerPart[100(t-3600 ore-60 minuti-secondi)];
	fm[ore]<>":"<>fm[minuti]<>":"<>fm[secondi]<>"."<>fm[decimi]
]
Moto[MR_,MC_,k_,v_,θ_,s0_]:=
Module[{M0=MR+MC,V,s,t0,eq,g=(9.8/1000),M,τx,τy,τ,t1,fig,x,y,tf,sol,ϵ=(10^(-6))},
M[t_]=M0-k t;
  (* Equazione del moto sulla rampa : 
M[t] s''[t] = k v - M[t] g Sin[θ] con {s[0]=0, s'[0]=0}; l'equazione può essere integrata direttamente senza usare DSolve *)
V[t_]=v Log[(M0/(M[ t]))]-g Sin[θ] t;     s[τ_]=Module[{t},Simplify[(∫_0^τ)(v Log[(M0/(M[t]))]-g Sin[θ] t)\[DifferentialD]t]];
(* calcoliamo l'istante di uscita dalla rampa *)
     t0=t/.FindRoot[s[t]==s0,{t,0.3}][[1]];
(* Condizioni iniziali per il moto successivo  *)
     x0=Re[s[t0]Cos[θ]];
     y0=Re[s[t0]Sin[θ]];
     vx0=Re[V[t0]Cos[θ]];
     vy0=Re[V[t0]Sin[θ]];
(* Istante in cui si esaurisce il carburante: *)
     t1=MC/k;
(* Equazioni del moto fino all'esaurimento del carburante: *)
τx[t_]:=((x'[t])/(Sqrt[((x'[t])^2)+((y'[t])^2)]));
τy[t_]:=((y'[t])/(Sqrt[((x'[t])^2)+((y'[t])^2)]));
eq={x''[t]==((k v τx[t])/(M[t])),y''[t]==((k v τy[t])/( M[t]))-g,x[t0]==x0,y[t0]==y0,x'[t0]==vx0,y'[t0]==vy0};
     sol=NDSolve[eq,{x,y},{t,t0,t1},Method->ExplicitRungeKutta,MaxSteps->10000][[1]];(* In conclusione le equazioni parametriche della traiettoria nei vari intervalli di tempo sono date da: *)
{x[t_],y[t_]}={x[t],y[t]}/.sol;
     x[t_]:=s[t]Cos[θ]/;t≤t0;
     y[t_]:=s[t]Sin[θ]/;t≤t0;(* Finito il carburante (t>t1) il razzo si muove soggetto alla sola forza peso *)
     x[t_]:=x[t1]+(t-t1) x'[t1]/;t>t1;
     y[t_]:=y[t1]+(t-t1) y'[t1]-(1/2)g (((t-t1))^2)/;t>t1;      tf=t1+((y'[t1]+(Sqrt[((y'[t1])^2)+2 g y[t1]]))/g);
  fig=ParametricPlot[{x[t],y[t]},{t,0,tf},DisplayFunction->Identity];
     Print["Istante del distacco dalla rampa = ",t0,
" s\nIstante (Subscript[t, 1]) in cui si esaurisce il carburante = ",SetPrecision[t1,3]," s\nx[(Subscript[t, 1])] = ",SetPrecision[x[t1],3],   " Km;  y[(Subscript[t, 1])] = ",SetPrecision[y[t1],3]," Km; \nx'[(Subscript[t, 1])] = ",SetPrecision[x'[t1],3]," Km/s;  y'[(Subscript[t, 1])] = ",SetPrecision[y'[t1],3]," Km/s\nGittata = ",SetPrecision[ x[tf],4]," Km"];    	
Table[
	Show[Graphics[Text[time[t],{0.8x[tf],5}]],Graphics[{AbsolutePointSize[7],Point[{x[t],y[t]}]}],ParametricPlot[{x[τ],y[τ]},{τ,0,t},PlotStyle->Dashing[{0.02,0.02}],DisplayFunction->Identity],PlotRange->{{-0.01,x[tf]+1},{-0.01,6}},Axes->True,AspectRatio->1,AxesLabel->{"x(Km)","y(Km)"}],
	{t,ϵ,tf+ϵ,tf/20}]
];

(Subscript[M, r])=700 ;
(Subscript[M, c]) =2100;
k=54.5;
θ=(π/4);
v=2.4 ;
t0=(((Subscript[M, c]) )/k);
d=0.01;
fig=Moto[(Subscript[M, r]),(Subscript[M, c]) ,k,v,θ,d]
fig1=ListAnimate[fig]

Export["N:\MyDir\Didattica\Fisica Computazionale\Mathematica 6.0\proiettile.gif",fig,"GIF"]

7.5.7 Problema: Pendolo.

Studiamo il moto di un corpo appeso a un filo lungo l. La forza esercitata dal filo è nulla se il filo non è teso; elastica, con costante elevata, se è teso: se (Overscript[(Subscript[x, 1]), ->])= {x1 , y1} sono le coordinate del corpo appeso all'estremo libero del filo (l'altro estremo è fisso nell'origine del sistema di coordinate), la forza del filo vale:
(Overscript[(Subscript[F, e l]), ->]) = - k (|(Overscript[(Subscript[x, 1]), ->])| - l ) ((Overscript[(Subscript[x, 1]), ->])/(|(Overscript[(Subscript[x, 1]), ->])|))       con    k > 0       se |(Overscript[(Subscript[x, 1]), ->])| > l
(Overscript[(Subscript[F, e l]), ->]) = 0      altrimenti.

Sul corpo agisce la forza (Overscript[(Subscript[F, d]), ->])  più la forza di gravità  m (Overscript[g, ->]) :

(Overscript[F, ->]) = (Overscript[(Subscript[F, e l]), ->]) +  m (Overscript[g, ->])

Clear[f,Pendolo,FilmPendolo,m,l,r0,v0,t0,t1]
Pendolo[f_,m_,l_,r0:{_,_},v0:{_,_},t0_,t1_]:=Pendolo[f,m,l,r0,v0,t0,t1]=
Module[{r,x,y,t},
r[t_]={x[t],y[t]};
eq=Flatten[Map[{#[[1,1]]==#[[2,1]],#[[1,2]]==#[[2,2]]}&,{m r''[t]==f[m,l,r[t],r'[t]],r[0]==r0,r'[0]==v0}]];
{x,y}={x,y}/.NDSolve[eq, {x,y}, {t, t0, t1},MaxStepSize->.1,StartingStepSize->.1,MaxSteps->100000][[1]]
];
FilmPendolo[t0_,t1_,np_,d_]:=
Module[{dt=((t1-t0)/np)},
Manipulate[
curva=ParametricPlot[{x[t],y[t]},{t,t0,t},DisplayFunction->Identity,PlotStyle->RGBColor[1,0,0]];
Show[
Graphics[{RGBColor[1,0,0],PointSize[.04],Point[{x[t],y[t]}]}],Graphics[Line[{{0,0},{x[t],y[t]}}]],
curva,Axes->True,PlotRange->d,AspectRatio->Automatic],
{t,t0+(10^(-6)),t1,dt}]
];

Dati e condizioni iniziali:

Clear[F,x,y];
F[m_,l_,r:{_,_},v:{_,_}]:=Module[{g=9.8`,K=(10^5),d=(Sqrt[r.r])},{If[(Sqrt[r.r])>l,-K (1-(l/(Sqrt[r.r]))) r[[1]],0],If[(Sqrt[r.r])>l,-K (1-(l/(Sqrt[r.r]))) r[[2]]-m g,-m g]}];
m=1.`;
l=10;
t0=0.`;
t1=30.`;
r0={0.365` l Sin[(π/4)],-l Cos[(π/4)]};
v0={0,0};
{x,y}=Pendolo[F,m,l,r0,v0,t0,t1];
Plot[x[t],{t,t0,t1},PlotRange->All]
Plot[y[t],{t,t0,t1},PlotRange->All]

np=200;
d={{-10,10},{0,-13}};
FilmPendolo[t0,t1,np,d]