Funktioner för att räkna med (ko-)vektorfält, system- transformationer samt att lösa system av första ordningens PDE:er

Erik Frisk

Version: 2004-03-04

BeginPackage["BasicGeo`"]

Diff::usage = "Compute the differential/gradient of a function \
λ w.r.t. a set of variables\n
Diff[λ,variables]"
LieDerivative::usage = "Compute the Lie-derivate of either a scalar \
λ along a vector field f, or a co-vector field ω along a \
vector field f\n
LieDerivative[λ,f,vars]\n
LieDerivative[ω,f,vars]"
JacobianMatrix::usage = "Compute the Jacobian of a vector field\n
JacobianMatrix[f,vars]"
LieProduct::usage = "Compute the Lie-product between two vector \
fields\n
LieProduct[f,g,vars]"
ChangeCoordinates::usage = "Change coordinates in a control-affine \
non-linear system.\n
x'=f(x)+g(x)u\n
y=h(x)\n
z=Φ(x)\n
ChangeCoordinates[f,g,h,Φ,vars]\n
\n
The vectorfield g and function h are optional."
Flow::usage = "Compute the flow of a vectorfield f, with independent \
variable t, starting at point x0\n
Flow[f,t,x0,vars]\n
Note: Highly experimental"
FrobeniusInverse::usage = "Compute the inverse transformation \
obtained when following the flows in F\n
FrobeniusInverse[F,x0,origvars,newvars]\n
Note: Highly experimental"
FrobSolve::usage = "Compute all solutions to the Frobenius PDE\n
FrobSolve[F,x0,vars,d]\n
Note: Highly experimental"
DistMemberQ::usage = "Determine if a vectorfield f is a member of a \
non-singular distribution Δ=span{f1,...,fd}. DistMemberQ[f,\
Δ]"

Begin["`Private`"]

Grundläggande funktioner

I detta avsnitt definiera funktionerna Diff[], LieDerivative[], JacobianMatrix[] samt LieProduct[], ChangeCoordinates[], DistMemberQ[].

Differential/gradient, Jakobian av ett vektorfält

Diff[λ_, vars_] := 
        Map[D[λ, #] &, vars];
Diff[f_List, vars_] := Map[Diff[#, vars] &, f];

JacobianMatrix[f_,  vars_] := Diff[f, vars];

LieDerivata

LieDerivative[λ_,  f_, vars_]  :=  
        Diff[λ,  vars] . f;
LieDerivative[w_List,  f_, 
          vars_]  :=  Simplify[{f} . 
                Transpose[JacobianMatrix[w, vars]] + 
{w} . 
                JacobianMatrix[f, vars]][[1]];

LieProduct

LieProduct[f_, g_,  vars_] := 
        JacobianMatrix[g,  vars] . f  -  
          JacobianMatrix[f,  vars] . g;

ChangeCoordinates

ChangeCoordinates[f_, Φ_, oldvar_] := 
        Module[{newsys}, 
newsys = 
            ChangeCoordinates[f, 0*f, 0, Φ, 
              oldvar]; 
Return[
            newsys[[1]]]];

    
ChangeCoordinates[f_, g_, Φ_, oldvar_] := 
        Module[{newsys}, 
newsys = 
            ChangeCoordinates[f, g, 0, Φ, 
              oldvar]; 
Return[
            newsys[[{1, 2}]]]];

    
ChangeCoordinates[f_, g_, h_, Φ_, oldvar_] := 
        Module[{J, Φinv, f2, g2, h2, 
            newvar}, 
newvar = Φ /. {Rule[x_, 
                    y_] -> 
                  x}; 

 (* 
            Compute Jacobian *) 
J = 
            JacobianMatrix[newvar /. Φ, 
              oldvar]; 

 (* 
            Invert transformation .  
                Assume unique inverse exists *) 
Φinv  = Solve[(# == (# /. Φ)) & /@ 
                  newvar, 
oldvar][[1]]; 

 (* 
            Convert f *) 
f2 = 
            Simplify[
              J . f /. Φinv]; 

 (* 
            Convert g .  Map over all g : 
                s if more than one *) 
g2 = 
            If[Head[g[[1]]] === List, 
Map[
                Simplify[J . # /. Φinv] &, 
                g], 
Simplify[
                J . g /. Φinv]]; 

 (* 
            Convert h .  Map over all h : 
                s if more than one *) 
h2 = 
            If[Head[h] === List, 
Map[
                Simplify[# /. Φinv] &, 
                h], 
Simplify[
                h /. Φinv]]; 

Return[{f2, g2, h2}];
];

DistMemberQ

DistMemberQ[f_, Δ_] := 
        Module[{n, d, z, rhs, eq}, 
d = 
            Length[Δ]; 
n = 
            Length[f]; 
consts = (C[#] & /@ 
                Range[1, d]); 
rhs = 
            consts . Δ; 
eq = (f[[#]] == rhs[[#]] & /@ Range[1, n]); 
Return[
            Length[Solve[eq, consts]] > 0]];

Frobeniusrelaterade funktioner

Flow

Flow[f_, z_, x_, 
          vars_] := 
DSolveWithX0[
            CreateFlowDiff[f, vars], 
x, 
            vars] /. {t -> z};

FindConstants[sol_, x0_] := 
        Module[{solconstants}, 
 (* 
            Hitta vilka konstantern från ODE - 
              lösaren som finns med *) 
solconstants = 
            Select[Variables[sol], 
              MatchQ[#, 
                  C[n_Integer]] &]; 

 (* Invertera nu systemet och skicka tillbaka resultatet *) 
Return[Solve[x0[[#]] == sol[[#]] & /@ 
                  Range[1, Length[x0]], 
                solconstants][[1]]];
];

CreateFlowDiff[f_, vars_] := 
        Module[{totvars, 
            primevars}, 
totvars = # -> #[t] & /@ 
              vars; 
primevars = #'[t] & /@ 
              vars; 
Return[primevars[[#]] == (f[[#]] /. totvars) & /@ 
              Range[1, Length[vars]]];
];

DSolveWithX0[diffeq_, x0_, vars_] := 
        Module[{totvars, diffsol, 
            consts}, 
totvars = # -> #[t] & /@ 
              vars; 
  
diffsol = (vars /. 
                  totvars) /. (DSolve[diffeq, vars /. totvars, 
                    t][[1]]); 
  
consts = 
            FindConstants[diffsol /. {t -> 0}, 
              x0]; 
Return[
            diffsol /. consts];
];

FrobeniusInverse

FrobeniusInverse[F_, x0_, origvars_, 
          newvars_] := 
Module[{x, 
            n}, 
n = 
            Length[origvars]; 
x = 
            x0; 
For[k = n, 
            k > 0, k--, 
x = 
                Flow[F[[k]], newvars[[k]], origvars, 
                    origvars] /. (origvars[[#]] -> 
                          x[[#]] & /@ 
                      Range[1, 
                        n]);
]; 
Return[Simplify[x]];
];

FrobSolve

Global`FrobSilent = True;

FrobSolve[F_, x0_, vars_] := 
        FrobSolve[F, x0, vars, 0];
FrobSolve[F_, x0_, vars_, d_] := 
        Module[{n, 
            foovars, Ψ, Ψinv}, 
If[
            Global`FrobSilent, 
            Off[Solve::incnst, Solve::ifun, 
              InverseFunction::ifun]]; 
n = 
            Length[F[[1]]]; 
foovars = z[#] & /@ 
              Range[1, n]; 
Ψ = 
            FrobeniusInverse[F, x0, vars, 
              foovars]; 
Ψinv = Solve[vars[[#]] == Ψ[[#]] & /@ Range[1, n], 
                foovars][[1]]; 
If[Global`FrobSilent, 
            On[Solve::incnst, Solve::ifun, 
              InverseFunction::ifun]]; 
Return[
            Take[foovars /. Ψinv, -(n - 
                    d)]];
];

Paket: Postamble

End[];

EndPackage[];