Scientific Calculator

(public) gio/AlgebraMatriciale

By gio gio

funzioni di algebra lineare

Tagged: matrici vettori

display=function(x,titolo){confirm(string(x),titolo)};

ID=identity(3);
// tensore antisimmetrico di Ricci Curbastro tridimensionale
Ricci=[[[0.,0.,0.],[0.,0.,1.],[0.,-1.,0.]],
	[[0.,0.,-1.],[0.,0.,0.],[1.,0.,0.]],
	[[0.,1.,0.],[-1.,0.,0.],[0.,0.,0.]]];

//"[k,j,flag,M(k,j)]=nozero(M[,err])"
// primo elemento != 0 ricercato per righe
nozero=function(M,err){
	var er=10;if err==null -> er=10 else er=err;
	var (k,j,f)=[0,0,true];
	while (f and (k< dim(M)[0]))->(k+=1;j=0;
        while (f and (j< dim(M)[1]))->(j+=1;f=ifzero(M(k,j),er)));
    [k,j,not f,M(k,j)]};
	
//danno come risposta il vettore (riga o colonna) di posizione "k"
estrairow=function(M,k){M[k-1]''};
estraicol=function(M,k){M'[k-1]'};
// eliminano una riga e colonna in una matrice (la numerazione comincia da 1)
eliminarow=function(v,j){var k=0;var d=dim(v)[0];var v1=[];
    while (k<=d-1 )-> (if (k != j-1) -> (v1= concat(v1,[v[k]]));k+=1);
	v1};
eliminacol=function(v,j){eliminarow(v',j)'};

//combinando le due funzioni superiori trova il minore M(j,l) di una matrice
// non cambia il segno secondo la parità degli indici
minor=function(v,j,l){eliminacol(eliminarow(v,j),l)};

//inseriscono la matrice "N" alla riga (o colonna) "k" della matrice "M".
// "k"=[1...dim(M)+1]
insertrow=function(M,N,k){var(j,d,f,O)=[0,dim(M)[0],false,[]];
	while (j<d)-> (
		if ((j != (k-1)) or f) -> (O=concat(O,[M[j]]);j+=1);
		if ((j) == (k-1))->(O=concat(O,N));f=true);
	O};
insertcol=function(M,N,k){insertrow(M',N',k)'};

//moltiplicano una riga (o colonna) per un numero
multrow=function(M,r,x){var N=clone(M);N[r-1]=N[r-1]*x;N};
multcol=function(M,c,x){var N=clone(M);N=N';N[c-1]=N[c-1]*x;N'};

//scambiano tra loro due righe (o due colonne)
scambiarow=function(M,k,j){var N=clone(M);
    var o=N[k-1];N[k-1]=N[j-1];N[j-1]=o;N};
scambiacol=function(M,k,j){var N=clone(M);N=N';
    var o=N[k-1];N[k-1]=N[j-1];N[j-1]=o;N'};

//calcola la inversa di una matrice,senza dividere per il determinante della matrice origine
//che potrebbe essere zero
pseudoinversa=function(M){var (k,j,d)=[0,0,0];
	if (dim(M)[0]==dim(M)[1]) -> (d=dim(M)[0]) 
	else (error("Matrice non quadrata;function:pseudoinversa"));
	var O=zeros(d,d);
	while (k<d)-> (while (j<d) -> ((O[j])[k]= det(minor(M,k+1,j+1))*(-1)^(j+k);j+=1);
	k+=1;j=0);
	O};

// da la matrice coniugata
CONJ=function(M){var(N,j,k)=[clone(M''),0,0];
	while (j< dim(N)[0])->
		(k=0;
		while (k<dim(N)[1])-> (N[j][k]=conj(N[j][k]);k+=1);
		j+=1;);
	N};

//prodotto vettoriale di due vettori tridimensionali
prvet=function(a,b){aa=clone(a);bb=clone(b);aa=aa';bb=bb';
	if (dim(aa)[1]==3)-> aa=aa';if (dim(bb)[1]==3) -> bb=bb';
	 pseudoinversa(augment([0;0;0],aa,bb))[0]'};
	 
//prodotto scalare di due vettori (matrici colonna)
prsc=function(a,b){aa=clone(a);bb=clone(b);aa=aa';bb=bb'';
	(aa*CONJ(bb))(1,1)};

//modulo e versore di un vettore
MOD=function(v) {real(sqrt((prsc(v,v))))};
VERS=function(v){if ifzero(MOD(v))-> error("vettore nullo;function;VERS") else v/MOD(v)};

//seno e coseno tra due matrici vettore
COS=function(v,w){var p=0.;p=prsc(v,w)/(MOD(v)*MOD(w));p};
SIN=function(a,b){var (m,n,f,ss)=[MOD(a),MOD(b),true,0];
	if ifzero(m*n)->(f=false;)
		else (ss=MOD(prvet(a,b))/(m*n));
	[ss,f]};

// calcola il volume del parallelepipedo formato dai versori di 3 vettori
// se zero i 3 vettori sono complanari
// se > zero formano una terna sinistrorsa, se < zero destrorsa
complanari=function(a,b,c){det(augment(VERS(a),VERS(b),VERS(c)))};

//trova il seno tra due vettori; se zero i due vettori sono allineati.
//in tal caso k è il rapporto con segno tra i due vettori
diplin=function(v1,v2){
	var (v,j,f,k)= [prvet(v1,v2),1,true,0];
	while (ifzero(v2(j,1)) and (j<=3)) ->(j+=1);
	if (j<=3) -> (k=v1(j,1)/v2(j,1));
	[|VER(v)|,k]};
// simile alla precedente
inlinea=function(a,b,c){diplin((b-a),(c-a));
	var (v,j,f,k)= [prvet(v1,v2),1,true,0];
	while (ifzero(v2(j,1)) and (j<=3)) ->(j+=1);
	if j>3-> f=false else(k=v1(j,1)/v2(j,1));
	[(f and ifzero(MOD(v))),k]};

// angolo orientato tra le proiezioni dei vettori "a" e "b" sul piano perpendicolare a "c"
tors=function(a,b,c){var (a,b)=[prvet(a,c),prvet(b,c)];var p=prvet(a,b);
	COS(c,p)*acos(COS(a,b))};

// costruisce la matrice di rotazione tridimensionale con asse "asse", angolo "angolo"
// la rotazione è anti-oraria se vista dalla punta del vettore "asse"
// "rd" assume i valori "d"=degrè, "g"= gradi centesimali , altro = radianti
Mrot=function(asse,angolo,rd){aa=clone(asse);ang=clone(angolo);
	if (rd=="d" or rd=="D")  -> (ang=ang*pi/180);
	if (rd=="g" or rd=="G")  -> (ang=ang*pi/200);
	if (system:type(aa)=="list") -> (aa=aa');
	aa=VERS(aa);
	var (x,y,z)=[aa(1,1),aa(2,1),aa(3,1)];
	var av=[e^(2i*pi/3);1;e^(-2i*pi/3)];
	if  (not ifzero(|x-y|+|y-z|)) ->
			av=[1-x^2-x*(y+z)-(y-z)*1i,
			1-y^2-y*(z+x)-(z-x)*1i,
			1-z^2-z*(x+y)-(x-y)*1i]';av=VERS(av);
	var U=augment(aa,av,CONJ(av));
	var AUTOV=Mdiagonal([1,e^(ang*1i),e^(-ang*1i)]);
	U*AUTOV*(1/U)};
	
//costruisce la matrice di autovettori colonna della matrice "M" e gli autovalori nella liista "a"
MMM=function(M,a){M*Mdiagonal(a)*(1/M)};

//rango di una matrice qualunque
RANGO=function(M,err){var er=10;if err==null -> er=10 else er=err;
		cerca=function(lst,r,err){var er=10;if err==null -> er=10 else er=err;
			var (ls,j,f)=[clone(lst),0,false];
			while (j<r) ->(shift(ls);j+=1;);
			map(x ->(if (ifzero(x,er) and (not f))->(j+=1) else (f=true)),ls);
			[j+1,f]};
	var (N,k,j,zr,fl)=[clone(M),0,0,0,true];var t=[0,false];
	if dim(M)[0]>dim(M)[1] ->N=N';var (r,c)=dim(N);
	while (k<r) -> (
		fl=true;
		if ifzero(N(k+1,k+1),er)->(t=cerca(N[k],k,er);
			if t[1] -> (N=scambiacol(N,k+1,t[0]))
			else (t=cerca(N'[k],k,er);
				if t[1] -> (N=scambiarow(N,k+1,t[0]))
				else (fl=false;zr+=1))
		);
		if fl -> (j=k+1;while (j<r) ->(N[j] += -N[k]*N[j][k]/N[k][k];j+=1;));
		k+=1);
	[r-zr,N]};
	
// costruisce la matrice diagonale da una lista "lst" di numeri
Mdiagonal=function(lst){var k=dim(lst)[0];var O=zeros(k,k);k=0;
    for j in lst -> (O[k][k]=j;k+=1);
    O};

// costruisce in una lista: la matrice degli autovettori normalizzati (per colonna),
// la matrice diagonale degli autovalori, la lista degli autovalori,la matrice "M"
// (non funziona per matrici non diagonalizzabili; solo per matrici 3X3)
AUTOvet=function(M,err){var er=10;if err==null -> er=10 else er=err;
	ATT=function(M){var (r,c,fl,va)=nozero(M,er);
		if fl -> VERS(M*conj(va))else
			error("Vettore Nullo;function: ATT in AUTOvett")};
    AT=function(M){var PS=pseudoinversa(M);
        var (r,c,fl,va)=nozero(PS,er);
        if fl -> VERS(estraicol(PS,c)*conj(va))else
			error("Matrice Nulla;function: AT in AUTOvett")};
	MA=function(aut,err){
		var au=clone(aut);var K=zeros(dim(M)[0],1);
		var ty =sort(tally(au,(x,y)->(ifzero(x-y,er))),(x,y)->(x[1]-y[1])<0 ? -1:1);
		if (dim(ty)[0]==1)-> (identity(dim(au)[0]))
		else (
		if (dim(ty)[0]== dim(au)[0])->(
			map(x-> (K=augment(K,AT(M-ID*x))),au);
			eliminacol(K,1))
		else (au=ty'[0];
			var x0=AT(M-ID*au[0]);
			var (x1,x2)=[zeros(dim(M)[0],1),zeros(dim(M)[0],1)];
			while ifzero(MOD(x1)*MOD(x2)*MOD(prvet(x1,x2)),er) ->(
				x1[0][0]=random([2])-1;x1[1][0]=random([2])-1;x1[2][0]=random([2])-1;
				x2[0][0]=random([2])-1;x2[1][0]=random([2])-1;x2[2][0]=random([2])-1;
				if ifzero(ty(1,2),er) -> (x1=x1-M*x1;x2=x2-M*x2) 
					else (x1=(M-ID*au[0])*x1;x2=(M-ID*au[0])*x2);
			x1=ATT(x1);x2=ATT(x2-prsc(x2,x1)*x1));
			autoval=[ty(1,1),ty(2,1),ty(2,1)];
			augment(x0,x1,x2)))};
	var tr=sum(diagonal(M));
	var mm=sum(diagonal(pseudoinversa(M)));
	var autoval=eq3([-det(M),mm,-tr,1]);
	var O=MA(autoval,er);
	var nn=Mdiagonal(autoval);
	[O,nn,autoval,M]};

// trova il punto di incontro tra la retta per i punti "r1,r2" e il piano di vettore direttore "d"
// e termine noto "tn" -> (<p,d>+tn=0 o altrimenti ( ax+by+cz+tn=0)
RP=function(r1,r2,d,tn){
	var k=(prsc(r1,d)+tn)/prsc(r1-r2,d);
	(1-k)*r1+k*r2};

// come sopra. il piano però è dato da tre suoi punti
RP1=function(r1,r2,p1,p2,p3){RP(r1,r2,prvet((p1-p3),(p2-p3)),-prsc(prvet(p2,p3),p1))};

// costruisce la matrice = f(M); la matrice "M" deve essere diagonalizzabile
//NOTA: non applica la funzione "f" a tutti gli elementi della matrice "M",
// vedi esempio sotto
fM=function(f,M,err){var(U,D,auto)=AUTOvet(M,err);
	U*Mdiagonal(map(f,auto))*(1/U)};

//due esempi applicativi: exp(M) e M^ex
// per la soluzione per sistemi lineari di equazioni differenziali del primo ordine
expM=function(M){
	expMM=function(x){e^x};
	fM(expMM,M)};
// esempio: la matrice N=potenzaM(H,1/3) è tale per cui H=N*N*N
potenzaM=function(M,ex,err){
    potenzaMM=function(ex,x){potenza(x,ex)};
    fM(bind(potenzaMM,ex),M,err)};

// date le rette p=k1*d1+r1 p=k2*d2+r1 trova la minima distanza e i valori corrispondenti [k1,k2]
// d1 e d2 sono i vettori direzione, r1,r2 i punti di passaggio
DMin=function(r1,d1,r2,d2){var (k1,k2)=[0,0];
		fk=function(ra,da,rb,db){
		var k=prsc(da*MOD(db)^2-db*prsc(db,da),ra-rb)/(((MOD(da)*MOD(db))^2)-prsc(db,da)^2)};
	if MOD(prvet(d1,d2))==0 ->  (k1=COS(r2-r1,d1))
	else [MOD(k1*d1+r1-k2*d2-r2),fk(r1,d1,r2,d2),fk(r2,d2,r1,d1)]};

spam? | offensive?

0 Comments

Sign in to leave a comment