Jacobian aus FindRoot extrahieren?

Question

Ich muss ein System von Festpunktgleichungen lösen und dann die Eigenwerte des Jacobi am Festpunkt berechnen. Es gibt ungefähr 50 Gleichungen mit 50 Variablen und sie enthalten viele numerische Integrale. Es wäre wirklich schwierig für mich, ein explizites Arbeitsbeispiel zu geben. Ich habe das Gleichungssystem in der Form

eqnSys={expr1,expr2,expr3,..,expr50};

Der Punkt ist, es ist eine Liste von Ausdrücken, die am Ende gleich Null sein muss und nicht als Funktion definiert ist, aber ich habe die Variablen in gespeichert

eqnVars={x1,x2,..,x50};

Ich habe eine sehr gute anfängliche Vermutung für die Wurzelfindung:

eqnGuess={{x1,1},{x2,2},..,{x50,50}};

Die grundlegende Newton-Methode von FindRoot bricht sofort zusammen und beschwert sich über das singuläre Jacobian. Also das

FindRoot[eqnSys,eqnGuess]

funktioniert nicht. Ich habe lange Zeit die Sekantenmethode verwendet, um die Wurzeln zu finden:

eqnGuessSec={{x1,0.9,1.1},{x2,1.8,2.2},..,{x50,45,55}};
FindRoot[eqnSys,eqnGuessSec]

Kürzlich bin ich auf die AffineCovariantNewton-Methode für FindRoot gestoßen, die wie ein Zauber wirkt und die Sekantenmethode zeitlich um den Faktor 4 übertrifft:

FindRoot[eqnSys,eqnGuess,Method -> {"AffineCovariantNewton"}]

Dem Bewertungsmonitor nach zu urteilen, hat es mehrere jakobianische Bewertungen:

FindRoot[eqnSys,eqnGuess,Method -> {"AffineCovariantNewton"},Jacobian -> {Automatic, EvaluationMonitor :> Print["J evaluated here"]}]

Meine Frage ist: Es wäre wirklich lukrativ für mich, den Jacobianer direkt vor FindRoot retten zu können. Ist es möglich, die von FindRoot erstellte Jacobi-Matrix zu extrahieren? Ich denke an so etwas

Reap@FindRoot[eqnSys,eqnGuess,Method -> {"AffineCovariantNewton"},Jacobian -> {Automatic,Sow[jacobian]}]

Ich interessiere mich nur für die rein numerische Matrix und nicht für eine symbolische. Bonusfrage: Was ist der effizienteste Weg, um das Gleichungssystem in eine Funktion umzuwandeln? Also so etwas wie

FeqnSys[x1_,x2_,...,x50_]:=eqnSys

Bearbeiten: Ich habe eine sehr einfache Version des Problems implementiert. Durch Erhöhen von UTrunc wird die Anzahl der Gleichungen erhöht (es werden dann jedoch zusätzliche Anfangsbedingungen benötigt). Ich brauche im Grunde das Objekt mit dem Namen ineedthisguy. Ich hoffte, dass es ohne diese analytische Unterscheidung erhalten werden kann, da ich für ein echtes Problem aufgrund der Speicherbeschränkung nur die vollständige Matrix in Blöcken erzeugen kann.

d = 3;
WorPrec = 16;
\[Alpha] = 1;

UTrunc = 6;

Z[r_] := 0
W[r_] := 0

U[r_] := Sum[
   ToExpression["u" <> ToString[n]]/n! (r - \[Kappa])^n, {n, 2, 
    UTrunc}];
\[Omega][r_] := U'[r] + 2 r U''[r]

MasterKernel1[d_, n1_, \[Omega]_?NumericQ, w_?NumericQ] := 
 MasterKernel1[d, n1, \[Omega], 
   w] = -2 \[Alpha] NIntegrate[
    E^-y y^(-1 + d/
      2) (1 + y) (y + w y^2 + E^-y \[Alpha] + \[Omega])^-n1, {y, 
     0, \[Infinity]}, 
    Method -> {Automatic, "SymbolicProcessing" -> False}, 
    WorkingPrecision -> WorPrec]

Derivative[1][MasterL[n_, d_]][\[Rho]_] := 
 Derivative[1][
   MasterL[n, 
    d]][\[Rho]] = -n (MasterL[n + 1, d][\[Rho]] \[Omega]'[\[Rho]] + 
     MasterL[pa][n + 1, d + 2][\[Rho]] Z'[\[Rho]] + 
     MasterL[pa][n + 1, d + 4][\[Rho]] W'[\[Rho]])

MasterL[n_, d_][\[Kappa]] := 
 MasterL[n, d][\[Kappa]] = 
  MasterKernel1[d, n, 2 \[Kappa] u2, W[\[Kappa]]]

BetaU[r_] := -d U[r] + (d - 2) r U'[r] - 
  1/(4 \[Pi]^2) MasterL[1, d][r]

dExpr[f_, betafunc_, n_] := D[k D[f[r], k] == betafunc[r], {r, n}]

GenBeta[f_, betafunc_, min_, max_] := Block[{expr, result, tmpres},
   expr = dExpr[f, betafunc, min];
   result = {(expr /. r -> \[Kappa])};
   Do[
    expr = D[expr, r];
    tmpres = Block[{r = \[Kappa]}, expr];
    result = Join[result, {tmpres}];
    , {i, min + 1, max}
    ];
   Return[result];
   ];

listU = GenBeta[U, BetaU, 1, UTrunc];
listU[[1]] = Thread[-listU[[1]]/u2, Equal];

FPEqn = ((Flatten@(List @@@ Flatten[listU]))[[2 ;; ;; 2]]);

varTrf = {g_[n_] :> ToExpression[ ToString[g] <> ToString[n]]};
varList = Flatten[{\[Kappa], Table[u[i], {i, 2, UTrunc}]}];

iniGuess = 
  Rationalize[
   List @@@ {\[Kappa] -> 0.04174875412610417566172053373396096686`12.,
      u2 -> 6.14584037490485804822706857376675685878`12., 
     u3 -> 60.04918116532118965443749174665446530096`12., 
     u4 -> 390.9010607033057646222`12., 
     u5 -> -3513.6112140902988423965`12., 
     u6 -> -93676.7079827356649900999`12.}, 0];
(*real solution:
{\[Kappa]\[Rule]0.0726928522670547`,u2\[Rule]4.570711765672155`,u3\
\[Rule]28.871831592476088`,u4\[Rule]134.9966784017132`,u5\[Rule]-371.\
15673934569224`,u6\[Rule]-14195.11815231752`}
*)

fOPT = Experimental`OptimizeExpression[FPEqn, 
   "OptimizationLevel" -> 2]; (*im not sure if this helps*)

fpLocator[initial__] := 
  FindRoot[fOPT // First, List @@@ initial, 
   Method -> {"AffineCovariantNewton"}, WorkingPrecision -> WorPrec, 
   StepMonitor :> {Print[initial[[All, 1]]]} ] ;

sol = fpLocator[iniGuess]


\!\(\*SuperscriptBox[\(MasterKernel1\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[d_, n1_, \[Omega]_?NumericQ, 
  w0_?NumericQ] := 
\!\(\*SuperscriptBox[\(MasterKernel1\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[d, n1, \[Omega], 
   w0] = -n1 MasterKernel1[d, n1 + 1, \[Omega], w0]

\!\(\*SuperscriptBox[\(MasterKernel1\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[d_, n1_, \[Omega]_?NumericQ, 
  w0_?NumericQ] := 
\!\(\*SuperscriptBox[\(MasterKernel1\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[d, n1, \[Omega], 
   w0] = -n1 MasterKernel1[d + 4, n1 + 1, \[Omega], w0]

ineedthisguy = Eigenvalues[D[FPEqn, {varList /. varTrf}] /. sol]

Answer 1

Hier ist ein Weg, um den Jacobian zu extrahieren. Die Idee ist, den linearen Löser in Sowden Jacobian und dann in Reapdiesen neu zu schreiben :

Richten Sie ein einfaches Problem ein:

f[X_] := Block[{x, y}, {x, y} = X; {Exp[x - 2] - y, y^2 - x}]
vars = {x, y};
start = {1, 1};
callFindRoot[f_, vars_, start_, opts___] := 
 vars /. FindRoot[f[vars], Evaluate[Transpose[{vars, start}]], opts]

Schreiben Sie das LinearSolveran Sowden Jakobianer neu:

MyLinearSolver = (Sow[#]; LinearSolve[##]) &;

Anruf FindRoot

Reap[callFindRoot[f, vars, start, Method -> {"AffineCovariantNewton"
    , "LinearSolver" -> {MyLinearSolver}
    (*,"BroydenUpdates"\[Rule]False*)
    }]]


(* {{0.019026, 
  0.137935}, {{{{0.367879, -1.}, {-1., 2.}}, {{0.295112, -1.}, {-1., 
     1.7796}}, {{0.191064, -1.}, {-1., 
     1.28831}}, {{0.096665, -1.}, {-1., 0.121763}}}}} *)

Ich werde dieses Beispiel zur Dokumentation hinzufügen.

Um die Gleichungen in eine Funktion umzuwandeln, können Sie Folgendes verwenden:

cf = With[{vars = vars, fun = f[vars]},
  Compile[{{X, _Real, 1}},
   Block[vars,
    vars = X;
    fun
    ]
   ]
  ]