Scientific Calculator

gio/AlgebraLineare

//NOTE: 07/03/2020
// I vettori sono inseriti come matrici colonna

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

// Crea il tensore di Ricci Curbastro di dimensione "k"	   
RICCI=function(k){
    IND=function(x,k){var (s,st)=[[],[]];var j=0; var cc= clone(x);var g=0;
        while (cc !=[]) -> (g=shift(cc);
            j=0;while (j<k) -> (push(st,(g mod k)); g= ipart(g/k);j+=1);
            s=push(s,reverse(st));st=[]);
            st=map(z-> (g=UZ(z);g==0)? 0 :(-1)^ordina((t,w)->(t<w),z)[0],s);st};
    var (j,c,s)=[0,[],[]];
    var ret=0;var h=-k;                 
    while (j<k) -> (push(c,j);j+=1);
    j=0;while (j<k)->(ret=MAP(ww->((j<(k-1)) ? c : (h+=k;IND(h+c,k))),ret);j+=1); 
[c,ret]};

//"[k,j,flag,M(k,j)]=nozero(M[,err])"
// primo elemento != 0 ricercato per righe
nozero=function(M,err){
	var er=(err==null ? 10 : 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 x
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};

// vettore di "superficie"
// la matrice M è costituita da n-1 vettori colonna di dimensione n
// Il vettore risultato è perpendicolare a tutti i vettori della matrice
// ATTENZIONE alla operazione di coniugazione: in R non fa differenza, ma 
// in C vi sono due metriche. Qui si è adottata la metrica ordinaria.
// Nella pseudo-metrica basta coniugare il risultato.
// vedi anche le due funzioni seguenti.
vetsup=function(M){CONJ(pseudoinversa(augment(zeros(dim(M)[0],1),M))[0]')};

//prodotto vettoriale di due vettori colonna tridimensionali
prvet=function(a,b){CONJ(pseudoinversa(augment([0;0;0],a,b)))[0]'};
	 
//prodotto scalare di due vettori (matrici colonna)
prsc=function(a,b){(a'*CONJ(b))(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){prsc(v,w)/(MOD(v)*MOD(w))};
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(aa,bb,cc){var (a,b)=[prvet(aa,cc),prvet(bb,cc)];var p=prvet(a,b);
	-COS(cc,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 stringa "d"/"D"/"g"/"G"/"r"/"R" secondo l' unità di misura dell' angolo
// se 'rd = null'  opp. altro -> si intende per default "radians"
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)]);
	CHOP(U*AUTOV*(1/U))};
	
//analoga alla precedente
MROT=function(asse,angolo,rd){
	var(c,s,xx,yy,zz,ang)=[0,0,clone(asse),[0;0;0]',[0;0;0]',clone(angolo)];
	if (rd=="d" or rd=="D") -> (ang=ang*pi/180);
	if (rd=="g" or rd=="G") -> (ang=ang*pi/200);  
	var s=e^(ang*1i);var c=(s+1/s)/2;s=(s-1/s)/2i;
    var RR=[1,0,0;0,c,-s;0,s,c];
	yy=prvet(xx,CONJ(xx))*1i;c=chop(MOD(yy));
	if (c==0) -> yy= VERS(CONJ(prvet(xx,[random()-1/2,random()-1/2,random()-1/2]')))
    else yy=VERS(yy);
	zz= VERS(CONJ(prvet(xx,yy)));
	var T=augment(xx,yy,zz);
	T*RR*T^-1};
	
//costruisce la matrice di autovettori colonna della matrice "A" e gli autovalori nella lista "a"
MMM=function(A,a){A*Mdiagonal(a)*(1/A)};

//costruisce la matrice simile a J tramite la matrice T
// Sia S=MT(J,T) -> S=T*J*T^-1
MT=function(J,T){T*J*(1/T)};

//rango di una matrice qualunque
RANGO=function(M,err){var er=(err==null ? 10 : err);
		cerca=function(lst,r,err){var er=(err==null ? 10 : 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};
	
// Forma di Jordan di una matrice 3x3 (anche numeri complessi)
// sia T=JORDAN(M) -> M=T[0]*T[1]*T[0]^-1
// costruisce in una lista:
//   T[0]: la matrice degli autovettori colonna(se possibile orto-normalizzati )
//   T[1]: la matrice diagonale di Jordan,
//   T[2]: la lista degli autovalori,
//   T[3]: restituisce la matrice "M"
JORDAN=function(M,err){
	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 JORDAN")};
	VRS=function(v,a,er){not ifzero(v(a,1),er) ? VERS(conj(v(a,1))*v):v};
	MA=function(aut,er){
			var au=clone(aut);var K=zeros(dim(M)[0],1);
			var (ra,x,y,z,d)=[0,[0,0,0]',[0,0,0]',[0,0,0]',[0,0,0]'];
			var ty =sort(tally(au,(x,y)->(ifzero(x-y,er/2))),(x,y)->(x[1]-y[1])<0 ? -1:1);
		TRIPLO=function(){var N=M-tr*ID/3;var ra=RANGO(N,er-2)[0];autoval=[tr,tr,tr]/3;
			rango1=function(){var r=nozero(N,er)[0];
					z=VRS(CONJ(N[r-1]'),3,er); y=N*z; x=VRS(CONJ(prvet(y,z)),1,er);augment(x,y,z)};
			rango2=function(){var r=nozero(N*N,er)[0];
					z=VRS(CONJ((N*N)[r-1]'),3,er);y=N*z;x=N*y;augment(x,y,z)};		
			CASE([ra==0,x->ID;ra== 2,rango2;true,rango1]);
			};	
		distinti=function(){var k=-1;
				map(x-> (K=augment(K,AT(M-ID*x))),au);K=eliminacol(K,1);
				map(x->(k+=1;ifzero(x[k],er) ? x : conj(x[k])*x/|x[k]|),K')'};
		doppio=function(){aa=ty(1,1);autoval=[aa,(tr-aa)/2,(tr-aa)/2];
				var N1=M-ID*aa;var N2=(M-ID*autoval[1])*N1*2/(tr-3*aa); ra=RANGO(N2,er-3)[0];	
				x=AT(N1);               			
				if (ra==1) ->(var r=nozero(N2,er)[0];z=VRS(N1*CONJ(N2[r-1]'),3,er);y=N2*z);								
				else
					(z=d;while ifzero(MOD(prvet(z,y),er))->(z=N1*MAP(x->(random()-1/2),d);y=N1*MAP(x->(random()-1/2),d));
					d=CONJ(prvet(y,z));
					if not ifzero(MOD(prvet(d,y),er)) -> (z=VRS(prsc(y,x)*z-prsc(z,x)*y,3,er));
					y=CONJ(VRS(prvet(d,z),2,er)));
				augment(x,y,z)};		
		CASE([dim(ty)[0]==1,TRIPLO;dim(ty)[0]== 3,distinti;true,doppio])};
	var er=(err==null ? 10 : err);		
	var tr=sum(diagonal(M));
	var autoval=eq3([-det(M),sum(diagonal(pseudoinversa(M))),-tr,1],er);
	autoval=sort(autoval,(x,y)-> (if (instring(system:type(CHOP(x)),"integer-float-rational")[0]) -> -1 else 1));
	var O=MA(autoval,er);
	[O,(1/O)*M*O,autoval,M]};

// esponenziale di matrice tramite sviluppo in serie di Taylor.
EXPM=function(M){var k=1;var S=identity(dim(M)[0]);var T=clone(S);delta=|MAXABS(M)|*10^-18;
    while (|MAXABS(T)|>delta)->(T=T*M/k;S=S+T;k+=1);
    S};
	
// 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

// per matrici diagonabilizzabili
fM=function(f,M,err){var(U,J,auto)=JORDAN(M,err);
	U*Mdiagonal(map(f,auto))*(1/U)};

// tre esempi applicativi: expM(M) LNM(M)e M^ex
// per la soluzione per sistemi lineari di equazioni differenziali del primo ordine
expM=function(M){fM(x->e^x,M)};
LNM=function(M){fM(x->ln(x),M)};
// esempio: la matrice N=potenzaM(H,1/3) è tale per cui H=N*N*N
potenzaM=function(M,ex){ fM(x->x^ex,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))};

// date le rette p=k1*d1+r1 p=k2*d2+r2 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 k=[0;0];
	var A = [MOD(d1)^2,-prsc(d1,d2);-prsc(d1,d2),MOD(d1)^2];
    var b= [prsc(d1,r2-r1);prsc(d2,r1-r2)];
	if det(A)==0 ->  [|SIN(r2-r1,d1)|/MOD(d1),[prsc(r2-r1,d1)/MOD(d1);0]]
	else (k=(1/A)*b;[MOD(k(1,1)*d1+r1-k(2,1)*d2-r2),k])};
spam? | offensive?
4 Comments
FirstTimeUser 10/13/2020 at 3:41 AM

What must one do to load this program, into what, and how to initiate it to run, and then enter a specific input to demonstrate that it is running as expected by the author?
Sonic the Hegehog 10/19/2023 at 5:45 PM

ummmmmmmmmmmmmm
???????????????????
Sonic the Hegehog 10/25/2023 at 2:49 PM

yo what the heke = is i sign not just = ok?
Sonic the Hegehog 10/25/2023 at 2:49 PM

yo what the heke = is i sign not just = ok?
gio/AlgebraLineare (function() { var po = document.createElement('script'); po.type = 'text/javascript'; po.async = true; po.src = 'https://apis.google.com/js/plusone.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(po, s); })();

gio/AlgebraLineare
