(*https://users.fit.cvut.cz/~stampfra/lectures/lin/lin-prednaska-14-handout.v1.pdf*)

(*popovidat o experimentalni matematice s pouzitim nahodnych cisel... slo by o tom delat doktorat na statistice a pravdepodobnosti... ale jistotu inzenyrum neda nic, tedy pravdepodobnost a statistika zvetsi pravdepodobnost, ze jsme nevypocitali pitomost :-]*)

Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
mf=MatrixForm;
sortByImag:=Sort[#,Im[#1]<Im[#2]&]&;
randomMatrix[rnm_]:=Table[RandomReal[{-1,1}]&[i],{i,rnm},{j,rnm}];

SeedRandom[4];
aijMax=1;
rnm=3;
A=randomMatrix[rnm];
%//mf;
{lambdas,vectors}=Eigensystem[A]//N;
P=Transpose[vectors];
diag=Inverse[P].A.P//Chop;
%//mf
{sortByImag[Eigenvalues[diag]]-sortByImag[lambdas]}//Chop

(*functions of matrices...*)
SeedRandom[4];
aijMax=1;
rnm=3;
A=randomMatrix[rnm];
fuctionOfMatrix[matrix_,function_]:=Module[{lambdas,vectors,P,invP,diagonalMatrix},
{lambdas,vectors}=Eigensystem[matrix]//N;
P=Transpose[vectors];
invP=Inverse[P];
diagonalMatrix=Inverse[P].matrix.P;
P.DiagonalMatrix[function/@Diagonal[diagonalMatrix]].invP
];
fuctionOfMatrix[A,#&]-A//Chop
fuctionOfMatrix[A,((#)^2)&]-A.A//Chop
fuctionOfMatrix[A,((#)^(-1))&]-Inverse[A]//Chop

(*dokazte, ze pokud takto zavedeme funkce, ktere maji konvergentni rozvoje v Taylorovy rady, tak ze vztahy platne pro funkce jedne realne (var: komplexni)plati i pro prislusne funkce "rozumnych" matic, ukazeme "dukaz" s pouzitim nahodne matice*)
sin[matrix_]:=fuctionOfMatrix[matrix,Sin];
cos[matrix_]:=fuctionOfMatrix[matrix,Cos];
exp[matrix_]:=fuctionOfMatrix[matrix,Exp];
ln[matrix_]:=fuctionOfMatrix[matrix,Log];
sqr[matrix_]:=fuctionOfMatrix[matrix,((#)^2)&];
rnm=3;
alfa=randomMatrix[rnm];
beta=randomMatrix[rnm];
mf/@{alfa,beta};

sin[2*alfa]-2*sin[alfa].cos[alfa]//Chop
exp[I*alfa]-(cos[alfa]+I*sin[alfa])//Chop
ln[exp[alfa]]-alfa//Chop
sqr[alfa]-exp[2*ln[alfa]]//Chop

tmax=10;
nPoints=40;
mat[t_]:=alfa*t;
vydejCislo[fce_,matSkalarFce_,t_]:=matSkalarFce@fuctionOfMatrix[mat[t],fce]//Re;
dejPlot[fce_,matSkalarFce_]:=Module[{dats},
dats={#,vydejCislo[fce,matSkalarFce,#]}&/@Range[0.,tmax,(tmax/nPoints)];
ListPlot[dats,AxesLabel->{"t","fce[t]"},Joined->True,PlotStyle->{Thick,Red},GridLines->Automatic]
];
matSkalarFce=(#/.List->Plus)&;
dejPlot[Sin,matSkalarFce]

(*co takhle vyzkouset, neni-li duvodem, ze nejsou stejne diagonalni matice a transformacni matice P... za body: najdete ruzne matice,ntere maji stejny Eigensystem a vyzkousejte identity pro dve promenne!*)
sqr[alfa+beta]-(sqr[alfa]+sqr[beta]+2*alfa.beta)//Chop
sqr[alfa+beta]-(sqr[alfa]+sqr[beta]+2*beta.alfa)//Chop
(*nepomuze ani prehodit poradi matic pri nasobeni!*)

(*pozor na zobecneni, pro vic promennych to nefunguje (alespon ne obecne)!!!*)
sin[alfa+beta]-(sin[alfa].cos[beta]+cos[alfa].sin[beta])//Chop;
sin[alfa+beta]-(sin[alfa].cos[beta]+sin[beta].cos[alfa])//Chop;
sin[alfa+beta]-(cos[beta].sin[alfa]+sin[beta].cos[alfa])//Chop;
sin[alfa+beta]-(cos[beta].sin[alfa]+cos[alfa].sin[beta])//Chop;

SeedRandom[1]
A=Table[RandomInteger[{-5,5}],{3},{3}]//N;
mf[A]
{lambdas,vectors}=Eigensystem[A];
P=Transpose[vectors];
invP=Inverse[P];
diag=invP.A.P;
mf[diag//Chop]
P.diag.invP//Chop//mf
sortByImag@Eigenvalues[diag]-sortByImag@lambdas//Chop

(*linear 1st order differential equations and eigenvalues*)
A=(({{2, 1}, {1, 2}}));
mrank=MatrixRank[A];
{lambdas,vectors}=Eigensystem[A];
P=Transpose[vectors];
invP=Inverse[P];
diag=invP.A.P;
mf@diag;
P.diag.invP//mf;
Eigenvalues[diag]-lambdas;
vars=(Subscript[x, #])[t]&/@Range[mrank];
eqns=Thread[D[vars,t]==A.vars];
vectInitConds=Range[mrank];
initConds=Thread[vars==vectInitConds]/.t->0;
expDiagMat=DiagonalMatrix[Exp[Diagonal[t*diag]]];
expDiagMat//mf;
expMatA=P.expDiagMat.invP;
sol=expMatA.vectInitConds//Simplify
DSolve[Union[eqns,initConds],vars,t]

(*"dukaz" ze rovnice n-teho radu ma stejne koreny jako matice vlastni cisla*)
ClearAll[a];
mrank=5;
origEqn=D[y[t],{t,mrank}]+(∑_(i=0)^(mrank-1))(Subscript[a, i])*D[y[t],{t,i}];
vars=(Subscript[x, #])[t]&/@Range[mrank];
eqns1=(Subscript[x, #-1])'[t]==(Subscript[x, #])[t]&/@Range[2,mrank];
eqnn=(Subscript[x, mrank])'[t]==-(∑_(i=1)^mrank)(Subscript[a, i-1])*(Subscript[x, i])[t];
eqnsTogether=Append[eqns1,eqnn]
{b,m}=CoefficientArrays[eqnsTogether,vars]//Normal
Thread[b==-m.vars]
A=-m;
mf@A
Det[A-λ*IdentityMatrix[mrank]]
origEqn/.{(Derivative[n_][y][t]):>(λ^n),y[t]->1}

(*matice a slozky*)
dosa={a->Exp[I*2*(Pi/3.)]};
T=(({{1, 1, 1}, {1, (a^2), a}, {1, a, (a^2)}}))/.dosa;
invT=(1/3)(({{1, 1, 1}, {1, a, (a^2)}, {1, (a^2), a}}))/.dosa;

Chop[Inverse[T]-invT//Simplify]

impedances={Zv,Z1,Z2};
matZ=RotateRight[impedances,#]&/@Range[0,2,1];
mf@matZ
transformedMatZ=Chop@Simplify[invT.matZ.T//Chop];
mf@%;

eqnJthPhase[j_]:=(Subscript[Δu, j])==(Subscript[R, j])*(Subscript[i, j])+(∑_(k=1)^3)(Subscript[Z, k,j])*(Subscript[i, k]);
eqns3phases=eqnJthPhase/@Range[3]//Simplify
varsToSolve=(Subscript[i, #])&/@Range[3];
{rhsabc,Zaux}=CoefficientArrays[eqns3phases,varsToSolve]//Normal;
Zabc=-1*Zaux;
mf@Zabc;
Thread[rhsabc==Zabc.varsToSolve]

symmetric={(Subscript[R, _]):>R,(Subscript[Z, k_,k_]):>Zv,(Subscript[Z, 1,2]):>Zm,(Subscript[Z, 3,2]):>Zm,(Subscript[Z, 1,3]):>Zm,(Subscript[Z, k_,r_]):>(Subscript[Z, r,k])};
ZabcSymmetric=Zabc//.symmetric;
mf@ZabcSymmetric;
Δvoltages=Thread[((Subscript[Δu, #])&/@Range[3])->ΔUf*{1,(a^2),a}];
values={ΔUf->100,Zm->20I,Zv->30I,R->1};
originalEqns=((eqns3phases//.symmetric)/.Δvoltages)/.values/.dosa
Solve[originalEqns][[1]]

Zcomponents=invT.(ZabcSymmetric/.values).T//Chop;
mf@Zcomponents;
voltagesOriginal=((Subscript[Δu, #])&/@Range[3]/.Δvoltages)/.values/.dosa;
voltagesComponents=invT.voltagesOriginal//Chop;
currentsComponents=(Subscript[iomn, #])&/@Range[3]
eqnsComponents=Thread[voltagesComponents==Zcomponents.currentsComponents]
solComponents=Solve[eqnsComponents][[1]]
originalCurrents=T.(currentsComponents/.solComponents)