> restart:
# HF, Übung 1: Lineare Algebra mit einem Computeralgebrasystem  (gewählt
# Maple)
# 
# Algebra der Vektoren und Tensoren, Spaltenmatrizen und quadratischen 
# Matrizen
# 
# Zuerst ist Einlesen des Pakets Lineare Algebra erforderlich:
# 
> with(linalg);
# 
# Addition
# 
> a:=vector(3);
> b:=vector(3):
> print(a);print(b);
> aplusb:=evalm(a+b);
# 
# Skalarprodukt
# 
> apunktb:=innerprod(a,b);
> norm(a,2);
> angle(a,b);
# 
# Vektorprodukt
# 
> akreuzb:=crossprod(a,b);
> c:=vector(3):
> bkreuzc:=crossprod(b,c);
# 
# Spatprodukt
# 
> spat1:=innerprod(akreuzb,c);
> 
> spat2:=innerprod(a,bkreuzc);
> simplify(spat1-spat2);
> matrixabc:=stackmatrix(a,b,c);
> spat3:=det(matrixabc);
> simplify(spat1-spat3);
> simplify(spat3-spat2);
# 
# Entwicklungssatz für zweifache Vektorprodukte
# 
> doppelt:=evalm(crossprod(a,bkreuzc)-innerprod(a,c)*b
> +innerprod(a,b)*c);
> map(simplify,doppelt);
# 
# Zahlenwerte
# 
> a:=vector([2,3,1]);
> b:=vector([-3,1,-2]);
> a[2];
> aplusb;
> print(aplusb);
> map(eval,aplusb);
> apunktb;
> map(eval,akreuzb);
> 
> restart:
# 
# Quadratische Matrizen
# 
> with(linalg):
> A:=matrix(3,3);
> B:=matrix(3,3);
> c:=vector(3):
> print(A); print(B);
> At:=transpose(A);
> evalm(A&*B);
> evalm(A&*c);
# 
# Invarianten
# 
> T:=matrix(3,3);
> trace(T);#Spur
> det(T);#Determinante
> id:=diag(1,1,1);#Einheitsmatrix
# 
# Die Eigenwertaufgabe
# 
> P:=evalm(T-sigma*id);
> dP:=det(P);
> pol:=sort(dP,sigma);
> polc:=collect(pol,sigma);
> c3:=coeff(polc,sigma,3);
> c2:=coeff(polc,sigma,2);
> c1:=coeff(polc,sigma,1);
> c0:=coeff(polc,sigma,0);
> c2-trace(T);
> c0-det(T);
# 
# Mit Zahlenwerten, symmetrischer Spannungstensor als Beispiel:
# 
> T:=matrix([[110.3,-22.8,13.1],[-22.8,-25.,-12.7],
> [13.1,-12.7,37.9]]);
# Ein Normaleneinheitsvektor:
> n:=normalize(vector([1.,1.,1.]));
# Der zugehörige Spannungsvektor:
> t:=evalm(T&*n);
# 
# Das Eigenwertproblem
# 
> P:=evalm(T-sigma*id);
> pol:=sort(det(P),sigma);
# Das ist die kubische Gleichung zur Bestimmung der drei Eigenwerte
> plot(pol,sigma=-50..150);
> ew:=fsolve(pol,sigma);
# 
# Wesentlich einfacher geht das in Maple so:
# 
> eww:=evalf(Eigenvals(T,Q));#Eigenwerte
> print(Q);#Spaltenmatrizen der Eigenvektoren
# 
# Kontrolle des Ergebnisses:
> for i to 3 do 
> evalm(T&*col(Q,i)-eww[i]*col(Q,i)) od;
# 
# Die Matrix  Q  ist orthogonal, also die Komponentenmatrix eines
# orthogonalen Tensors. Probe:
# 
> evalm(transpose(Q)&*Q);
> det(Q);
> iQ:=inverse(Q);
> tQ:=transpose(Q);
> evalm(iQ-tQ);
>