> restart:
# HF, Übung2:  Vektor- und Tensoranalysis
# 
# Operationen in kartesischen und Zylinderkoordinaten
# 
> with(linalg):
# 
# Bereitstellung zweier Prozeduren zur Berechnung von Vektorgradient und
# Tensordivergenz:
# 
> vecgrad:=proc()
>  local i,v,q,co,h,vq,hq,hinv,dia,vdia,ddia;
>  description ` Vektorgradient:  v dyadisch nabla`;
>  v:=args[1];
>  q:=args[2];
>  h:=[1,1,1];
>  if nargs>2 then 
>    co:=args[3];
>    if type(co,`=`)then
>      if rhs(co)=cylindrical then
>        h:=[1,q[1],1] fi;
>      if rhs(co)=spherical then
>        h:=[1,q[1],q[1]*sin(q[2])] fi
>    else h:=co
>    fi
>  fi;
>  vdia:=array(sparse,1..3,1..3):
>  for i to 3 do vdia[i,i]:=v[i] od;
>  vq:=jacobian(v,q);
>  hq:=jacobian(h,q);
>  hinv:=array(sparse,1..3,1..3):
>  for i to 3 do hinv[i,i]:=1/h[i] od;
>  dia:=evalm(hq&*hinv&*v);
>  ddia:=array(sparse,1..3,1..3):
>  for i to 3 do ddia[i,i]:=dia[i] od;
>  evalm((vq-hinv&*transpose(hq)&*vdia
>  +ddia)&*hinv);
>  end: 
> 
> tensdiv:=proc()
>  local i,l,t,q,co,h,hinv,hq,vol,tz,tzeile,di;
>  description ` Tensordivergenz:  t punkt nabla`;
>  t:=args[1];
>  q:=args[2];
>  h:=[1,1,1];
>  if nargs>2 then 
>    co:=args[3];
>    if type(co,`=`)then
>      if rhs(co)=cylindrical then
>        h:=[1,q[1],1] fi;
>      if rhs(co)=spherical then
>        h:=[1,q[1],q[1]*sin(q[2])] fi
>    else h:=co
>    fi
>  fi;
>  hinv:=array(sparse,1..3,1..3):
>  for i to 3 do hinv[i,i]:=1/h[i] od;
>  vol:=h[1]*h[2]*h[3];
>  tz:=evalm(t&*hinv*vol);
>  hq:=jacobian(h,q);
>  for i to 3 do
>      di[i]:=diverge(row(tz,i),q)/vol +1/h[i]*
>  add((t[l,i]*hq[i,l]-t[l,l]*hq[l,i])/h[l],
>  l=1..3)
>  od;
>  vector([di[1],di[2],di[3]])
>  end:
>  
# Abkürzungen:
> alias(zyl=([r,phi,z], coords=cylindrical)):
> alias(kart=([x,y,z])):
# 
# Definition eines skalaren Feldes  theta  (z.B. Temperatur) und 
# eines Vektorfeldes  u  (z.B. Verschiebung)
# in Abhängigkeit von kartesischen bzw. Zylinderkoordinaten:
# 
> Theta_k:=theta(x,y,z);
> Theta_z:=theta(r,phi,z);
> u_k:=vector([ux(x,y,z),uy(x,y,z),uz(x,y,z)]);
> u_z:=vector([ur(r,phi,z),uphi(r,phi,z),uz(r,phi,z)]);
# 
# Gradient des skalaren Feldes:
# 
> grad_Theta_k:=grad(Theta_k,kart);
> grad_Theta_z:=grad(Theta_z,zyl);
# 
# Überprüfung der Identität   rot grad = 0:
# 
> rotgrad_k :=curl(grad_Theta_k,kart);
> rotgrad_z :=curl(grad_Theta_z,zyl);
# 
# Anwendung von Delta = div grad  auf das skalare Feld
# 
> delta_Theta_k:=diverge(grad_Theta_k,kart);
> laplacian(Theta_k,kart);# Zum Vergleich
> delta_Theta_z:=expand(diverge(grad_Theta_z,zyl));
> expand(laplacian(Theta_z,zyl));# Zum Vergleich
# 
# Divergenz und Gradient des Vektorfeldes:
# 
> div_u_k:=diverge(u_k,kart);
> div_u_z:=expand(diverge(u_z,zyl));
> grad_u_k:=vecgrad(u_k,kart);
> grad_u_z:=vecgrad(u_z,zyl);
# 
# Probe, dass die Divergenz die Spur des Gradienten ist:
# 
> trace(grad_u_k)-div_u_k;
> simplify(trace(grad_u_z)-div_u_z);
# 
# Rotation des Vektorfeldes:
# 
> rot_u_k:=curl(u_k,kart);
# 
# Probe, dass die Rotation durch den antimetrischen Teil des Gradienten
# bestimmt ist (und umgekehrt):
# 
> grad_u_k[3,2]-grad_u_k[2,3]-rot_u_k[1];
# 
# Überprüfung der Identität   div rot = 0
# 
> diverge(rot_u_k,kart);
# 
# Dasselbe in Zylinderkoordinaten:
# 
> rot_u_z:=curl(u_z,zyl);
> simplify(grad_u_z[3,2]-grad_u_z[2,3]-rot_u_z[1]);
> diverge(rot_u_z,zyl);
# 
# Gradient der Divergenz des Vektorfeldes
# 
> graddiv_u_k:=grad(div_u_k,kart);
> graddiv_u_z:=grad(div_u_z,zyl);
# 
# Die Divergenz eines Tensorfeldes
# 
> T_k:=matrix([[txx(x,y,z),txy(x,y,z),txz(x,y,z)],
> [tyx(x,y,z),tyy(x,y,z),tyz(x,y,z)],
> [tzx(x,y,z),tzy(x,y,z),tzz(x,y,z)]]);
> 
> tensdiv(T_k,kart);
> T_z:=matrix([[trr(r,phi,z),trphi(r,phi,z),trz(r,phi,z)],
> [tphir(r,phi,z),tphiphi(r,phi,z),tphiz(r,phi,z)],
> [tzr(r,phi,z),tzphi(r,phi,z),tzz(r,phi,z)]]);
> 
> tensdiv(T_z,zyl);
# 
# Anwendung von Delta = div grad  auf das Vektorfeld
# 
> delta_u_k:=tensdiv(grad_u_k,kart);
> delta_u_z:=tensdiv(grad_u_z,zyl);
# 
# Zweimalige Anwendung der Rotation auf das Vektorfeld
# 
> rotrot_u_k:=curl(rot_u_k,kart);
> rotrot_u_z:=curl(rot_u_z,zyl);
# 
# Beweis, dass gilt   Delta = div grad    =   grad div  - rot rot :
# 
> evalm(delta_u_k-graddiv_u_k+rotrot_u_k);
> duz:=evalm(delta_u_z-graddiv_u_z+rotrot_u_z);
> map(simplify,duz);
# 
# Der Sonderfall rotationssymmetrischer Felder.
# Die Behandlung geschieht natürlich mittels Zylinderkoordinaten.
# 
> theta:=(r,phi,z)->g(r);
> map(eval,grad_Theta_z);
> delta_Theta_z;
> ur:=(r,phi,z)->uu(r);
> uphi:=0;
> uz:=0;
> map(eval,grad_u_z);
> map(eval,rot_u_z);
> div_u_z;
> map(eval,delta_u_z);
> map(eval,graddiv_u_z);
> 
>