<<JavaGraphics`

(* oscillatore smorzato e forzato *)

eqpendolo = {y''[t] + gamma y'[t] + omega0^2 y[t] == f0 Cos[omega t], y'[0] == 1, y[0]==0};


smorzato[gamma_, omega0_] = y[t] /. DSolve[ eqpendolo /. f0->0, y[t], t ];

forzato[f0_, omega_, omega0_]  = y[t] /. DSolve[ eqpendolo /. gamma->0, y[t], t ];

smorzatoforzato[gamma_, f0_, omega_, omega0_] = y[t] /. DSolve[ eqpendolo, y[t], t ];



grasmo1 = Plot[ Evaluate[ smorzato[0.1,1] ], {t,0,20}];
grasmo2 = Plot[ Evaluate[ smorzato[30,1] ], {t,0,20}];
smorzatocritico[omega0_] := Limit[ smorzato[gamma, omega0], gamma->2 omega0];
grasmo3 = Plot[ Evaluate[ smorzatocritico[1] ], {t,0,20}];
grasmorzato = Show[ grasmo1, grasmo2, grasmo3];


grafor1 = Plot[ Evaluate[ forzato[1,0.1,1] ], {t,0,100}];
grafor2 = Plot[ Evaluate[ forzato[1,10,1] ], {t,0,100}];
forzatocritico[f0_, omega0_] := Limit[ forzato[f0,omega, omega0], omega->omega0];
grafor3 = Plot[ Evaluate[ forzatocritico[1,1] ], {t,0,100}];
graforzato = Show[ grafor1, grafor2, grafor3];

grasmorzatoforzato = Plot3D[ Evaluate[ smorzatoforzato[1,1,omega,1] ], {t,0,100}, {omega,0.1,10}];




(******************
calcolo autovalori autovettori  nativo
 *)




A = {{-5, -9, 8}, {9, 3, -9}, {10, 2, 7}};
A = {{-5, -9, -9}, {9, 3, 8}, {10, 2, 7}};

eig = Eigenvalues[ A];
eigvec = Eigenvectors[ A];


(******************
calcolo autovalori autovettori  implementazione propria
 *)

caratteristico = Det[ A - lambda IdentityMatrix[3] ];
autovals = lambda /. Solve[ caratteristico==0, lambda];

componenti = Table[v[i], {i,1,3}];

autovecs = Table[ componenti /. Solve[ A.componenti == autovals[[j]] componenti, {v[1],v[2]}] /. v[3]->1 ,
		 {j,1,3}
		  ] // N;

(******************
calcolo autovalori autovettori  implementazione propria  caso generale NxN
 *)

  myEigvec[A_List] := Module[{lunghezza, caratteristico, autovals, componenti, autovecs, tmp},

			     lunghezza = Length[A];
			     caratteristico = Det[ A - lambda IdentityMatrix[lunghezza] ];
                             autovals = lambda /. Solve[ caratteristico==0, lambda];

                             componenti = Table[v[i], {i,1,lunghezza}];


                             autovecs = Table[ sol= Solve[ A.componenti == autovals[[j]] componenti, componenti];
					       tmp = componenti /. sol;
					       tmp,
		                              {j,1,lunghezza}
		                             ];
			     variabile = Union[ Cases[ autovecs, v[_], Infinity]][[1]];
			     autovecs = autovecs /. variabile ->1;
			     autovecs//N
			     
			     ];

(**** Gram-Schmidt ***)
base[vec_List] := Table[  
                            u[j] = vec[[j]] -
                                     Sum[ ( vec[[j]].u[k-1] )/( u[k-1].u[k-1] ) u[k-1] ,  {k,2,j}];


                            u[j] = u[j] / Sqrt[ u[j].u[j] ];

			    ris = u[j];
			    
                            ris,


                         {j,1,Length[vec]}
				     ];



ATA = Transpose[A].A;
V = N[ Eigenvectors[ ATA ] , 8 ];
Vort = base[ V ];
ATAeig = N[ Eigenvalues[ ATA ], 8 ];


S2 = N[ Table[ If[ i===j, ATAeig[[i]], 0 ], {i,1,3}, {j,1,3} ], 8 ];

S = Sqrt[S2];

check1 = ATA.V === V.S2;

U = A.Vort.Inverse[S];

check2 = A === U.S.Transpose[Vort];

(****************)


DerivaMonomio[expr_, var_] := Module[{risultato},

                                      If[  MatchQ[ expr, w_. var^n_.],
                                           prefattore = expr /. w_. var^n_. -> w;
                                           esponente = expr /. w_. var^n_. -> n;
                                           ];
                                      
                                      Which[ MatchQ[ expr, w_. var^n_.] && FreeQ[prefattore,var] &&
                                                  FreeQ[ esponente, var],  
                                                  risultato = prefattore esponente var^(esponente-1),

                                             FreeQ[ expr, var],
                                             risultato = 0,

                                             True,
                                             risultato = nonsoderivare[ expr, var]

                                            ];
                                      
                                      risultato
];



myDer[expr_, var_] :=   Module[ {polinomio, risultato},

                                polinomio = Expand[ expr];
                                Which[ Head[ polinomio] === Plus,
                                       risultato= Map[ DerivaMonomio[#,var ]&, polinomio] ,
                                       
                                       True,
                                       risultato = DerivaMonomio[ expr, var ]

                                       ];

                                risultato
                                
                                ];










 
