Estrarre Jacobian da FindRoot?

Question

Devo risolvere un sistema di equazioni di punto fisso e poi calcolare gli autovalori dello Jacobiano al punto fisso. Ci sono circa 50 equazioni con 50 variabili e contengono molti integrali numerici, sarebbe davvero difficile per me fornire un esempio operativo esplicito. Ho il sistema di equazioni nella forma

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

il punto è che è una lista di espressioni che alla fine deve essere uguale a zero e non è definita come funzione, ma ho le variabili salvate in

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

Ho un'ipotesi iniziale molto buona per la ricerca della radice:

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

Il metodo di base di Newton di FindRoot si rompe immediatamente, lamentandosi del singolare Jacobiano. Così questo

FindRoot[eqnSys,eqnGuess]

non funziona. Ho usato a lungo il metodo delle secanti per trovare le radici:

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

Recentemente mi sono imbattuto nel metodo AffineCovariantNewton per FindRoot che funziona come un fascino e supera il metodo secante nel tempo di un fattore 4:

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

A giudicare dal monitor di valutazione, ha diverse valutazioni giacobiane:

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

La mia domanda è: sarebbe davvero redditizio per me essere in grado di salvare il Jacobian direttamente da FindRoot. È possibile estrarre la matrice Jacobiana costruita da FindRoot? Sto pensando a qualcosa di simile

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

Mi interessa solo la matrice puramente numerica e non simbolica. Domanda bonus: qual è il modo più efficiente per trasformare il sistema di equazioni in una funzione? Quindi qualcosa di simile

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

Modifica: ho implementato una versione molto semplice del problema. Aumentando UTrunc si aumenta il numero di equazioni (ma sono necessarie ulteriori condizioni iniziali). Fondamentalmente ho bisogno dell'oggetto con il nome ineedthisguy. Speravo che potesse essere ottenuto senza questa differenziazione analitica, perché per un problema reale posso solo generare la matrice completa in blocchi a causa della limitazione della memoria.

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

Ecco un modo per estrarre lo Jacobiano. L'idea è di riscrivere il risolutore lineare Sowsullo Jacobiano e poi su Reapquello:

Imposta un semplice problema:

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]

Riscrivi il LinearSolveral Sowgiacobiano:

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

Chiamata 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}}}}} *)

Aggiungerò questo esempio alla documentazione.

Per trasformare le equazioni in una funzione potresti usare qualcosa sulla falsariga di:

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