Benutzer:Alva2004/Programme/VektorWxm

/* vektor algebra und analysis */
display2d: false$
grind=true$

load("eigen"); /* laed columnvector */

o:zeromatrix( 3,1 )$
O:zeromatrix( 3,3 )$
E:ident(3)$
ex:columnvector( [ 1,0,0 ] )$
ey:columnvector( [ 0,1,0 ] )$
ez:columnvector( [ 0,0,1 ] )$

/* =============== kartesische koord =============== */
spur(T):=mat_trace(T)$

postfix("^T")$
"^T" (MAT):=transpose(MAT)$

infix("dotT");
"dotT" (MAT1,MAT2):=mat_trace(MAT1^T.MAT2)$

/* kreuzprodukte */
infix( "xV" );
"xV" ( VEC1,VEC2 ):= columnvector(
    [ VEC1[ 2,1 ]*VEC2[ 3,1 ] - VEC1[ 3,1 ]*VEC2[ 2,1 ],
      VEC1[ 3,1 ]*VEC2[ 1,1 ] - VEC1[ 1,1 ]*VEC2[ 3,1 ],
      VEC1[ 1,1 ]*VEC2[ 2,1 ] - VEC1[ 2,1 ]*VEC2[ 1,1 ] ] )$

infix( "Tx" );
"Tx" ( MAT,VEC ):= ( ex.( ( MAT[ 1 ]^T xV VEC )^T )
                     + ey.( ( MAT[ 2 ]^T xV VEC )^T )
                     + ez.( ( MAT[ 3 ]^T xV VEC )^T ) )$

/* vi*ei x Tjk*ej o ek =  */
infix( "xT" );
"xT" ( VEC,MAT ):= ( ( VEC xV ( MAT^T[ 1 ]^T ) ).ex^T
                     + ( VEC xV ( MAT^T[ 2 ]^T ) ).ey^T
                     + ( VEC xV ( MAT^T[ 3 ]^T ) ).ez^T )$

/* auesseres tensorprodukt S # T nach formeltalg12.py, siehe cartbase.wxm */
infix( "xx" );
"xx" ( S,T ):= matrix([
S[2,2]*T[3,3]-S[2,3]*T[3,2]-S[3,2]*T[2,3]+S[3,3]*T[2,2],
-S[2,1]*T[3,3]+S[2,3]*T[3,1]+S[3,1]*T[2,3]-S[3,3]*T[2,1],
+S[2,1]*T[3,2]-S[2,2]*T[3,1]-S[3,1]*T[2,2]+S[3,2]*T[2,1] ],
[ -S[1,2]*T[3,3]+S[1,3]*T[3,2]+S[3,2]*T[1,3]-S[3,3]*T[1,2],
+S[1,1]*T[3,3]-S[1,3]*T[3,1]-S[3,1]*T[1,3]+S[3,3]*T[1,1],
-S[1,1]*T[3,2]+S[1,2]*T[3,1]+S[3,1]*T[1,2]-S[3,2]*T[1,1] ],
[ +S[1,2]*T[2,3]-S[1,3]*T[2,2]-S[2,2]*T[1,3]+S[2,3]*T[1,2],
-S[1,1]*T[2,3]+S[1,3]*T[2,1]+S[2,1]*T[1,3]-S[2,3]*T[1,1],
+S[1,1]*T[2,2]-S[1,2]*T[2,1]-S[2,1]*T[1,2]+S[2,2]*T[1,1] ]
)$

norm(v):=mat_norm(v,frobenius)$

I1(T):=mat_trace(T)$
I2(T):=1/2*( mat_trace(T)^2-mat_trace(T.T) )$
I3(T):=determinant(T)$

i( MAT ):=matrix( [ MAT[ 2,3 ] - MAT[ 3,2 ] ],
                  [ MAT[ 3,1 ] - MAT[ 1,3 ] ],
                  [ MAT[ 1,2 ] - MAT[ 2,1 ] ] )$

axialerTensor(v):=block(
    k:list_matrix_entries(v),
    x:k[1], y:k[2], z:k[3],
    matrix([0,-z,y],[z,0,-x],[-y,x,0] )
    );

drehTensor(v):=block(
    w:norm(v),
    ud:axialerTensor(unitvector(v)),
    ident(3) + sin(w)*ud + ( 1 - cos(w) )*ud.ud
    );

st:sin( t )$
ct:cos( t )$
tt:tan( t )$
sf:sin( f )$
cf:cos( f )$

/* =============== zylinderkoord =============== */
eR:columnvector( [ cf, sf, 0 ] )$
ef:columnvector( [ -sf, cf, 0 ] )$
zylV( VEC ):= trigsimp( columnvector( [ eR^T.VEC, ef^T.VEC, ez^T.VEC ] ) );
zylT( MAT ):= trigsimp( [ 
[ eR^T.MAT.eR, eR^T.MAT.ef, eR^T.MAT.ez ],
[ ef^T.MAT.eR, ef^T.MAT.ef, ef^T.MAT.ez ],
[ ez^T.MAT.eR, ez^T.MAT.ef, ez^T.MAT.ez ] ] );
divVzyl( VEC ):= trigsimp( eR^T.diff( VEC,R ) + 1/R*ef^T.diff( VEC,f )
                           + ez^T.diff( VEC, z ) )$
nablaDotTzyl( MAT ):= divVzyl( MAT )^T;
divTzyl( MAT ) := nablaDotTzyl( MAT^T )$
gradFzyl( ska ):= trigsimp( diff( ska,R )*eR + 1/R*diff( ska,f )*ef
                            + diff( ska,z )*ez )$
gradVzyl( VEC ):= trigsimp( diff( VEC,R ).eR^T + 1/R*diff( VEC,f ).ef^T
                            + diff( VEC,z ).ez^T )$
rotVzyl( VEC ):= trigsimp( eR xV diff( VEC,R ) + 1/R*ef xV diff( VEC,f )
                           + ez xV diff( VEC, z ) )$
nablaCrossTzyl( MAT ):= ( eR xT diff( MAT,R ) + 1/R*ef xT diff( MAT,f )
                          + ez xT diff( MAT,z ) )$
rotTzyl( MAT ):= nablaCrossTzyl( MAT^T )$
deltaFzyl( ska ):= ( diff( ska,R )/R + diff( ska,R,2 ) + diff( ska,f,2 )/R^2
                 + diff( ska,z,2 ) )$

/* =============== kugelkoordinaten =============== */
er:[ st*cf, st*sf, ct ]^T$
et:[ ct*cf, ct*sf, -st ]^T$
ef:[ -sf, cf, 0 ]^T$
sphV( VEC ):= trigsimp( columnvector( [ er^T.VEC, et^T.VEC, ef^T.VEC ] ) )$
sphT( MAT ):= trigsimp( [ 
[ er^T.MAT.er, er^T.MAT.et, er^T.MAT.ef ],
[ et^T.MAT.er, et^T.MAT.et, et^T.MAT.ef ],
[ ef^T.MAT.er, ef^T.MAT.et, ef^T.MAT.ef ] ] );
nablaDotVsph( VEC ):= trigsimp( er^T.diff( VEC,r )
                                + 1/r*et^T.diff( VEC,t )
                                + 1/(r*st)*ef^T.diff( VEC,f ) )$
divVsph( VEC ):=nablaDotVsph( VEC )$
gradFsph( SKA ):= ( trigsimp( diff( SKA,r )*er
                              + diff( SKA,t )/r*et
                              + diff( SKA,f )/( r*st )*ef ) )$
gradVsph( VEC ):= trigsimp( diff( VEC,r ).er^T
                            + 1/r*diff( VEC,t ).et^T
                            + 1/(r*st)*diff( VEC,f ).ef^T )$
rotVsph( VEC ):= trigsimp( er xV diff( VEC,r )
                           + 1/r*et xV diff( VEC,t )
                           + 1/(r*st)*ef xV diff( VEC,f ) )$

/* =============== kartesische koordinaten =============== */
/* divergenz eines vektorfeldes -> skalar vk*ek.Aij*eioej=Aij*vi*ej */
divV( VEC ):= ratsimp( ex^T.diff( VEC,x ) + ey^T.diff( VEC,y )
                     + ez^T.diff( VEC,z ) )$
nablaDotT( MAT ):= divV( MAT )^T; /* spaltenvektor */

divT( MAT ):=ratsimp( diff( MAT,x ).ex + diff( MAT,y ).ey + diff( MAT, z ).ez )$

/* gradient eines skalarfeldes -> spaltenvektor */
gradF( FKT ):= ratsimp( diff( FKT,x )*ex + diff( FKT,y )*ey + diff( FKT,z )*ez )$

/* gradient eines vektorfeldes -> tensor (N o v)^T */
gradV( VEC ):= ratsimp( diff( VEC,x ).ex^T + diff( VEC,y ).ey^T
                      + diff( VEC,z ).ez^T )$

/* rotation */
rotV( VEC ):= ex xV diff( VEC,x ) + ey xV diff( VEC,y ) + ez xV diff( VEC,z )$
nablaCrossT( MAT ):= ( ex xT diff( MAT,x ) + ey xT diff( MAT,y )
                     + ez xT diff( MAT,z ) )$
rotT( MAT ):= ( ex xT diff( MAT^T,x ) + ey xT diff( MAT^T,y )
              + ez xT diff( MAT^T,z ) )$

/* laplace */
deltaF( SKA ):= divV( gradF( SKA ) )$
deltaV( VEC ):= divT( gradV( VEC ) )$

ok:[];
T:genmatrix( T_, 3,3 )$
expand(charpoly(T,x))$
/*
lam:eigenvalues(T)$
evc:eigenvectors(T)$
*/