Ekstrak Jacobian dari FindRoot?

Question

Saya harus menyelesaikan sistem persamaan titik tetap dan kemudian menghitung nilai eigen dari Jacobian pada titik tetap. Ada sekitar 50 persamaan dengan 50 variabel dan mereka mengandung banyak integral numerik, akan sangat sulit bagi saya untuk memberikan contoh kerja eksplisit. Saya memiliki sistem persamaan dalam bentuk

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

intinya adalah, ini adalah daftar ekspresi yang pada akhirnya harus sama dengan nol dan tidak didefinisikan sebagai fungsi, tetapi saya memiliki variabel yang disimpan di

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

Saya memiliki tebakan awal yang sangat bagus untuk temuan root:

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

Metode Newton dasar FindRoot segera rusak, mengeluh tentang Jacobian tunggal. Jadi ini

FindRoot[eqnSys,eqnGuess]

tidak bekerja. Saya menggunakan metode garis potong untuk waktu yang lama untuk menemukan akarnya:

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

Baru-baru ini saya menemukan metode AffineCovariantNewton untuk FindRoot yang bekerja seperti pesona dan mengungguli metode garis potong dalam waktu dengan faktor 4:

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

Dilihat dari evaluasi monitor, ada beberapa evaluasi Jacobian:

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

Pertanyaan saya adalah: Akan sangat menguntungkan bagi saya untuk dapat menyelamatkan Jacobian langsung dari FindRoot. Apakah mungkin untuk mengekstrak matriks Jacobian yang dibangun oleh FindRoot? Saya sedang memikirkan sesuatu seperti

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

Saya hanya tertarik pada matriks numerik murni dan bukan matriks simbolik. Pertanyaan bonus: Apa cara paling efisien untuk mengubah sistem persamaan menjadi sebuah fungsi? Jadi seperti itu

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

Sunting: Saya telah menerapkan versi masalah yang sangat sederhana. Meningkatnya UTrunc bertambah sebagai banyaknya persamaan (tapi kemudian membutuhkan tambahan kondisi awal). Saya pada dasarnya membutuhkan objek dengan nama ineedthisguy. Saya berharap itu dapat diperoleh tanpa perbedaan analitis ini, karena untuk masalah nyata saya hanya dapat menghasilkan matriks penuh dalam potongan karena keterbatasan memori.

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

Berikut adalah cara untuk mengekstrak Jacobian. Idenya adalah untuk menulis ulang pemecah linier ke SowJacobian dan kemudian ke Reap:

Siapkan masalah sederhana:

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]

Tulis ulang LinearSolverke SowJacobian:

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

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

Saya akan menambahkan contoh ini ke dokumentasi.

Untuk mengubah persamaan menjadi fungsi, Anda dapat menggunakan sesuatu di sepanjang baris:

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