Extrair Jacobian de FindRoot?

Question

Tenho que resolver um sistema de equações de ponto fixo e então calcular os autovalores do Jacobiano no ponto fixo. Existem cerca de 50 equações com 50 variáveis ​​e elas contêm muitos integrais numéricos, seria muito difícil para mim dar um exemplo de trabalho explícito. Eu tenho o sistema de equações na forma

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

a questão é que é uma lista de expressões que tem que ser igual a zero no final e não é definida como uma função, mas eu tenho as variáveis ​​salvas em

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

Tenho uma estimativa inicial muito boa para a descoberta da raiz:

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

O método básico de Newton de FindRoot quebra imediatamente, reclamando sobre o singular Jacobiano. Então, é isso

FindRoot[eqnSys,eqnGuess]

não funciona. Usei o método da secante por muito tempo para encontrar as raízes:

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

Recentemente, descobri o método AffineCovariantNewton para FindRoot que funciona perfeitamente e supera o método secante no tempo por um fator de 4:

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

A julgar pelo monitor de avaliação, tem várias avaliações Jacobianas:

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

Minha pergunta é: seria muito lucrativo para mim ser capaz de salvar o Jacobiano diretamente de FindRoot. É possível extrair a matriz Jacobiana construída por FindRoot? Estou pensando em algo como

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

Estou interessado apenas na matriz puramente numérica e não simbólica. Pergunta bônus: Qual é a maneira mais eficiente de transformar o sistema de equações em uma função? Então algo como

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

Edit: implementei uma versão muito simples do problema. Aumentar UTrunc aumenta conforme o número de equações (mas precisa de condições iniciais adicionais). Basicamente, preciso do objeto com o nome ineedthisguy. Eu esperava que pudesse ser obtido sem essa diferenciação analítica, porque para um problema real eu só posso gerar a matriz completa em pedaços devido à limitação de memória.

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

Aqui está uma maneira de extrair o Jacobiano. A ideia é reescrever o solucionador linear para Sowo Jacobiano e depois para Reapele:

Configure um problema simples:

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]

Reescreva o LinearSolverpara Sowo Jacobiano:

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

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

Vou adicionar este exemplo à documentação.

Para transformar as equações em uma função, você poderia usar algo ao longo das linhas de:

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