Я пытаюсь в wolfram решить систему уравнений (24),...

python
import os, json, textwrap, re, pathlib, sys, subprocess, glob
paths = glob.glob('/mnt/data/*')
paths[:20], len(paths)

python
import os, glob, pathlib, re, textwrap, json, sys
paths=glob.glob('/mnt/data/*')
print(len(paths))
for p in paths:
    print(repr(p), os.path.getsize(p))

python
nb_path='/mnt/data/Гаусс-Бонне_new.nb'
data=open(nb_path,'rb').read()
print(data[:100])
text=data.decode('utf-8', errors='replace')
print(text[:1000])
print('len', len(text))

python
import shutil, subprocess, os, json, textwrap, re, sys, glob
print(shutil.which('wolframscript'), shutil.which('math'), shutil.which('Mathematica'))

python
for fp in ['/mnt/data/GB_shooting_fixed_v2.wl','/mnt/data/GB_shooting_fixed.wl']:
    print("FILE",fp)
    print(open(fp, encoding='utf-8').read()[:10000])

python
import numpy as np, math
from scipy.optimize import root_scalar
from scipy.integrate import solve_ivp

ell=0.2; mu=0.1; h=0.1; v=mu/h; rH=1.0; eps=1e-6; r0=1+eps; rmax=50
def V(x): return h*h/4*(x*x-v*v)**2
def Vp(x): return h*h*(x*x-v*v)*x
def b1(phiS,dphiS): return (1-rH*rH*V(phiS))/(rH-8*ell*ell*phiS*dphiS)
def horizonEq(phiS,dphiS):
    bb=b1(phiS,dphiS)
    return dphiS - 16*ell**2*phiS/rH**3 - 8*bb*ell**2*phiS/rH**2*(1+16*ell**2*phiS*dphiS/rH) - (Vp(phiS)+16*ell**2*phiS*V(phiS)/rH**2)/bb
def phiPrimeS(phiS):
    # scan roots?
    x0=24*ell**2*phiS/rH**3
    # use secant
    from scipy.optimize import root
    sol=root(lambda y: horizonEq(phiS,float(y[0])), [x0])
    if not sol.success:
        # try bracket
        xs=np.linspace(-10,10,1000)
        vals=[horizonEq(phiS,x) if abs(rH-8*ell**2*phiS*x)>1e-8 else np.nan for x in xs]
        for a,b,fa,fb in zip(xs[:-1],xs[1:],vals[:-1],vals[1:]):
            if np.isfinite(fa) and np.isfinite(fb) and fa*fb<=0:
                return root_scalar(lambda x:horizonEq(phiS,x), bracket=[a,b]).root
        raise Exception('no root')
    return sol.x[0]
def Acoef(r,pp,z,ll):
    return (1/r)*(1-(8*ell**2/r)*(1-3*np.exp(-ll))*pp*z)
def nuPrimeFun(r,pp,z,ll):
    return ((np.exp(ll)-1)/r**2 + z*z/2 - np.exp(ll)*V(pp))/Acoef(r,pp,z,ll)
def lambdaPrimeFun(r,pp,z,ll):
    return (-(np.exp(ll)-1)/r**2 - (16*ell**2/r**2)*(1-np.exp(-ll))*pp*z + z*z/2 + np.exp(ll)*V(pp))/Acoef(r,pp,z,ll)
def C33(r,pp,z,ll):
    q=np.exp(-ll); nup=nuPrimeFun(r,pp,z,ll)
    return 1 - (128*ell**4*pp**2*nup*q*(1-q))/(r**2*(r+16*ell**2*q*pp*z))
def rhs(r,y):
    pp,z,nn,ll=y
    q=np.exp(-ll); el=np.exp(ll)
    nup=nuPrimeFun(r,pp,z,ll)
    lamp=lambdaPrimeFun(r,pp,z,ll)
    a11=1 + (16*ell**2*q/r)*pp*z
    a12=(16*ell**2*q/r)*nup*pp
    btheta=(nup/2+1/r)*(nup-lamp) - (8*ell**2*q/r)*(3*nup*lamp - nup**2)*pp*z + (16*ell**2*q/r)*nup*z**2 + z**2 + 2*el*V(pp)
    a21=(8*ell**2*pp/r**2)*(1-q)
    a22=1
    bscalar=((nup-lamp)/2+2/r)*z + (8*ell**2*pp/r**2)*((1-q)*((nup-lamp)*nup/2)+q*nup*lamp) - el*Vp(pp)
    M=np.array([[a11,a12],[a21,a22]],float)
    nupp,phidd=np.linalg.solve(M, -np.array([btheta,bscalar]))
    return [z,phidd,nup,lamp]
def shoot(phiS, stop=True):
    d=phiPrimeS(phiS); bb=b1(phiS,d)
    if not np.isfinite(bb) or bb<=0: return None
    y0=[phiS+d*eps,d,math.log(eps),-math.log(bb*eps)]
    events=[]
    def evA(r,y): return Acoef(r,y[0],y[1],y[3])
    evA.terminal=True
    evA.direction=0
    def evC(r,y): return C33(r,y[0],y[1],y[3])
    evC.terminal=True
    def evTop(r,y): return y[0]-2*v
    evTop.terminal=True
    evTop.direction=0
    def evLow(r,y): return y[0]+v
    evLow.terminal=True
    evLow.direction=0
    events=[evA,evC,evTop,evLow]
    sol=solve_ivp(rhs, (r0,rmax), y0, method='Radau', events=events, rtol=1e-8, atol=1e-10, max_step=(rmax-r0)/1000)
    return sol
# test d
for ps in [0.001,0.1,0.5,0.6068,0.9]:
    try: 
        d=phiPrimeS(ps); print(ps,d,b1(ps,d),horizonEq(ps,d))
    except Exception as e: print('err',ps,e)

python
def residual_py(ps, allow_incomplete=True, verbose=False):
    try:
        sol=shoot(ps)
    except Exception as e:
        if verbose: print("fail",ps,e)
        return np.nan
    if sol is None:
        return np.nan
    val=sol.y[0,-1]-v
    if verbose:
        print("ps",ps,"tend",sol.t[-1],"phiend",sol.y[0,-1],"res",val,"status",sol.status, [arr for arr in sol.t_events if len(arr)])
    return val

for ps in np.linspace(0.01,0.99,10):
    print(ps, residual_py(ps, verbose=True))

python
from scipy.optimize import brentq
# Get bracket
xs=np.linspace(1e-8,0.999,80)
vals=[]
for x in xs:
    vals.append(residual_py(x))
for i in range(len(xs)-1):
    if np.isfinite(vals[i]) and np.isfinite(vals[i+1]) and vals[i]*vals[i+1]<=0:
        print("br",xs[i],vals[i],xs[i+1],vals[i+1])
        break
root=brentq(lambda x: residual_py(x), 0.55,0.67, xtol=1e-8)
root, residual_py(root, verbose=True)

python
import re, os, textwrap, json, pathlib
text=open('/mnt/data/Гаусс-Бонне_new.nb', encoding='utf-8', errors='replace').read()
# find occurrences
for pat in ['findPhiS', 'findBracket', 'residual', 'shoot', 'First@Select', 'br = First']:
    print(pat, text.find(pat), text.count(pat))
# show around findPhiS in raw
idx=text.find('findPhiS')
print(text[idx-2000:idx+3000])

python
for idx in [text.find('residual'), text.find('shoot')]:
    print('--- idx', idx)
    print(text[idx-1500:idx+2500])

python
idx=text.find('Cell[BoxData[\n RowBox[{"\\[IndentingNewLine]"')
# search around shoot more
idx=text.find('SHOOTING', 60000)
print(idx)
print(text[idx:idx+6000])

python
for pat in ['rS', 'eps', 'ell', 'mu', 'h', 'v']:
    print(pat, text.find(f'"{pat}"'))
# Show beginning code roughly
print(text[1000:8000])

python
# show around phiDDfun
for pat in ['phiDDfun', 'nuPrimeExpr', 'lambdaPrimeExpr']:
    idx=text.find(pat)
    print('---',pat,idx)
    print(text[idx-1200:idx+2500])

python
code = r'''(* ::Package:: *)

(* Fixed robust shooting for Babichev-Emond-Ramazanov,
   arXiv:2207.03944v2, Eqs. (24),(25),(27) + BC (3),(28),(32).

   How to run:
     1) Evaluation -> Quit Kernel -> Local
     2) Get["/full/path/GB_shooting_fixed_v3.wl"]

   The old notebook failed because:
     - rS was shifted to 1+eps instead of keeping the horizon at rH=1;
     - lambda was initialized as -Log[eps] instead of -Log[b1 eps];
     - residual returned artificial 10^6 for failed shots;
     - findPhiS used First[{}] when no sign-changing bracket was found;
     - phi'' was obtained from a differentiated divided expression.
*)

ClearAll["Global`*"];

wp = 50;
eps = SetPrecision[10^-6, wp];

rH = SetPrecision[1, wp];       (* true horizon radius *)
r0 = rH + eps;                  (* integration starts slightly outside *)
rmax = SetPrecision[50, wp];

ell = SetPrecision[0.2, wp];
mu  = SetPrecision[0.1, wp];
h   = SetPrecision[0.1, wp];
v   = mu/h;

V[x_?NumericQ]  := (h^2/4) (x^2 - v^2)^2;
Vp[x_?NumericQ] := h^2 (x^2 - v^2) x;

(* Eq. (30): e^(-lambda)=b1 (r-rH)+... *)
b1[phiS_?NumericQ, dphiS_?NumericQ] :=
  (1 - rH^2 V[phiS])/(rH - 8 ell^2 phiS dphiS);

(* Eq. (32), with the printed r evaluated at rH. *)
horizonEq[phiS_?NumericQ, dphiS_?NumericQ] := Module[{bb = b1[phiS, dphiS]},
  If[!NumericQ[bb] || bb == 0, Return[Indeterminate]];
  dphiS
   - (16 ell^2 phiS)/rH^3
   - (8 bb ell^2 phiS/rH^2) (1 + (16 ell^2 phiS dphiS)/rH)
   - (Vp[phiS] + (16 ell^2 phiS V[phiS])/rH^2)/bb
];

phiPrimeS[phiS_?NumericQ] := Module[{x, x0, fr},
  x0 = SetPrecision[24 ell^2 phiS/rH^3, wp];
  fr = Quiet@FindRoot[
      horizonEq[phiS, x] == 0,
      {x, x0},
      WorkingPrecision -> wp,
      AccuracyGoal -> 28,
      PrecisionGoal -> 28,
      MaxIterations -> 100
    ];
  If[Head[fr] =!= List, Return[$Failed]];
  x /. fr
];

Acoef[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, ll_?NumericQ] :=
  (1/rr) (1 - (8 ell^2/rr) (1 - 3 Exp[-ll]) pp zz);

nuPrimeFun[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, ll_?NumericQ] :=
  ((Exp[ll] - 1)/rr^2 + zz^2/2 - Exp[ll] V[pp])/
    Acoef[rr, pp, zz, ll];

lambdaPrimeFun[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, ll_?NumericQ] :=
  (-(Exp[ll] - 1)/rr^2
    - (16 ell^2/rr^2) (1 - Exp[-ll]) pp zz
    + zz^2/2 + Exp[ll] V[pp])/
    Acoef[rr, pp, zz, ll];

(* Eq. (33): zero of this coefficient is a genuine canonical singularity. *)
C33[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, ll_?NumericQ] := Module[
  {q = Exp[-ll], nup},
  nup = nuPrimeFun[rr, pp, zz, ll];
  1 - (128 ell^4 pp^2 nup q (1 - q))/(rr^2 (rr + 16 ell^2 q pp zz))
];

(* Use Eq. (26) and Eq. (27) as a 2x2 linear system for {nu'', phi''}.
   This is much more stable than differentiating nuPrimeFun. *)
rhsVals[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, nn_?NumericQ, ll_?NumericQ] :=
 Module[{q, el, nup, lamp, a11, a12, a21, a22, btheta, bscalar, sol2},
  q = Exp[-ll];
  el = Exp[ll];

  nup = nuPrimeFun[rr, pp, zz, ll];
  lamp = lambdaPrimeFun[rr, pp, zz, ll];

  If[!VectorQ[{q, el, nup, lamp}, NumericQ], Return[{Indeterminate, Indeterminate, Indeterminate, Indeterminate}]];

  (* Eq. (26): a11 nu'' + a12 phi'' + btheta == 0 *)
  a11 = 1 + (16 ell^2 q/rr) pp zz;
  a12 = (16 ell^2 q/rr) nup pp;
  btheta = (nup/2 + 1/rr) (nup - lamp)
    - (8 ell^2 q/rr) (3 nup lamp - nup^2) pp zz
    + (16 ell^2 q/rr) nup zz^2
    + zz^2 + 2 el V[pp];

  (* Eq. (27): a21 nu'' + a22 phi'' + bscalar == 0 *)
  a21 = (8 ell^2 pp/rr^2) (1 - q);
  a22 = 1;
  bscalar = ((nup - lamp)/2 + 2/rr) zz
    + (8 ell^2 pp/rr^2) ((1 - q) ((nup - lamp) nup/2) + q nup lamp)
    - el Vp[pp];

  sol2 = Quiet@Check[
     LinearSolve[{{a11, a12}, {a21, a22}}, -{btheta, bscalar}],
     {Indeterminate, Indeterminate}
   ];

  {zz, sol2[[2]], nup, lamp}
];

validSolutionQ[sol_] :=
  ListQ[sol] && Length[sol] == 4 &&
   And @@ (Head[#] === InterpolatingFunction & /@ sol);

shoot[phiS_?NumericQ, requireFull_: False] := Module[
  {pS, dphiS, bb, sol, dom},

  pS = SetPrecision[phiS, wp];

  dphiS = phiPrimeS[pS];
  If[dphiS === $Failed || !NumericQ[dphiS], Return[$Failed]];

  bb = b1[pS, dphiS];
  If[!NumericQ[bb] || bb <= 0, Return[$Failed]];

  sol = Quiet@NDSolveValue[
     {
      phi'[r] == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[1]],
      F'[r]   == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[2]],
      nu'[r]  == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[3]],
      lam'[r] == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[4]],

      phi[r0] == pS + dphiS eps,
      F[r0]   == dphiS,

      (* a1 is a time-normalization gauge; a1=1 is fine. *)
      nu[r0]  == Log[eps],

      (* Correct horizon expansion from Eq. (29),(30). *)
      lam[r0] == -Log[bb eps],

      WhenEvent[Acoef[r, phi[r], F[r], lam[r]] == 0, "StopIntegration"],
      WhenEvent[C33[r, phi[r], F[r], lam[r]] == 0, "StopIntegration"],

      (* Stop wrong shooting branches before they generate nDSolve::ndsz. *)
      WhenEvent[phi[r] - 2 v == 0, "StopIntegration"],
      WhenEvent[phi[r] + v == 0, "StopIntegration"],
      WhenEvent[Abs[F[r]] - 10 v == 0, "StopIntegration"]
     },
     {phi, F, nu, lam}, {r, r0, rmax},
     WorkingPrecision -> wp,
     AccuracyGoal -> 26,
     PrecisionGoal -> 26,
     MaxStepFraction -> 1/2000,
     Method -> {"StiffnessSwitching"}
   ];

  If[!validSolutionQ[sol], Return[$Failed]];

  dom = sol[[1]]["Domain"][[1, 2]];
  If[requireFull && dom < rmax - 100 eps, Return[$Failed]];

  sol
];

(* Residual:
   For a bad shot it is acceptable to use the last finite point of the
   InterpolatingFunction. That preserves the sign needed for bracketing. *)
residual[phiS_?NumericQ] := Module[{sol, dom, val},
  sol = shoot[phiS, False];
  If[sol === $Failed, Return[Indeterminate]];

  dom = sol[[1]]["Domain"][[1, 2]];
  val = Quiet@Check[sol[[1]][dom] - v, Indeterminate];

  If[NumericQ[val] && Im[N[val]] == 0, Re[val], Indeterminate]
];

findBracket[lo_: 10^-8, hi_: 0.999, n_: 250] := Module[
  {grid, vals, good, pairs, br},

  grid = SetPrecision[Subdivide[lo v, hi v, n], wp];

  vals = ({#, Quiet@Check[residual[#], Indeterminate]} & /@ grid);
  good = Select[vals, NumericQ[#[[2]]] && Im[N[#[[2]]]] == 0 &];

  If[Length[good] < 2,
    Print["No usable shooting samples. Try increasing wp, lowering rmax, or pass a manual bracket to findPhiS[{phi1,phi2}]."];
    Print["Raw sample count = ", Length[vals], ", usable = ", Length[good]];
    Return[$Failed]
  ];

  pairs = Partition[good, 2, 1];
  br = Select[pairs, #[[1, 2]] #[[2, 2]] <= 0 &];

  If[br === {},
    Print["No sign-changing bracket was found."];
    Print["Best sampled points {phiS,residual}: ",
      N[Take[SortBy[good, Abs[#[[2]]] &], UpTo[10]], 12]];
    Print["For the default parameters try explicitly: findPhiS[{0.55 v, 0.67 v}]"];
    Return[$Failed]
  ];

  First[br][[All, 1]]
];

findPhiS[br0_: Automatic] := Module[{br, x1, x2, r1, r2, fr},
  br = If[br0 === Automatic, findBracket[], SetPrecision[br0, wp]];
  If[br === $Failed, Return[$Failed]];

  {x1, x2} = br;
  r1 = Quiet@Check[residual[x1], Indeterminate];
  r2 = Quiet@Check[residual[x2], Indeterminate];

  If[!(NumericQ[r1] && NumericQ[r2]),
    Print["Manual bracket endpoints are not usable. {x1,r1,x2,r2} = ", N[{x1, r1, x2, r2}, 12]];
    Return[$Failed]
  ];

  If[r1 r2 > 0,
    Print["The bracket does not change sign. {x1,r1,x2,r2} = ", N[{x1, r1, x2, r2}, 12]];
    Print["Use scanResidual[] to inspect the shooting function."];
    Return[$Failed]
  ];

  fr = Quiet@FindRoot[
      residual[phiS] == 0,
      {phiS, x1, x2},
      WorkingPrecision -> wp,
      AccuracyGoal -> 20,
      PrecisionGoal -> 20,
      MaxIterations -> 100
    ];

  If[Head[fr] =!= List,
    Print["FindRoot failed. Last bracket data: ", N[{x1, r1, x2, r2}, 12]];
    Return[$Failed]
  ];

  phiS /. fr
];

scanResidual[n_: 30] := Module[{grid, vals},
  grid = SetPrecision[Subdivide[10^-8 v, 0.999 v, n], wp];
  vals = ({#, Quiet@Check[residual[#], Indeterminate]} & /@ grid);
  Grid[Prepend[N[vals, 12], {"phiS", "residual"}], Frame -> All]
];

(* For the default parameters this manual bracket avoids any ambiguity.
   If parameters are changed, either adjust it or use findPhiS[] to scan. *)
phiSsol = findPhiS[{SetPrecision[0.55 v, wp], SetPrecision[0.67 v, wp]}];

If[phiSsol === $Failed,
  Print["No solution was found by the current shooting setup."];
  Print["Run scanResidual[40] and choose two phiS values with opposite residual signs."];
  Abort[]
];

dphiSsol = phiPrimeS[phiSsol];
solFinal = shoot[phiSsol, True];

If[solFinal === $Failed,
  Print["Root was found, but final integration did not reach rmax."];
  Print["Try increasing wp or lowering rmax. phiS = ", N[phiSsol, 16]];
  Abort[]
];

Print["phiS = ", N[phiSsol, 16]];
Print["phiPrimeS = ", N[dphiSsol, 16]];
Print["phi(rmax)-v = ", N[solFinal[[1]][rmax] - v, 16]];
Print["domain end = ", N[solFinal[[1]]["Domain"][[1, 2]], 16]];

Print["min Acoef on grid = ",
  N[Min[Acoef[#, solFinal[[1]][#], solFinal[[2]][#], solFinal[[4]][#]] & /@
    Subdivide[r0, rmax, 1000]], 16]];

Print["min C33 on grid = ",
  N[Min[C33[#, solFinal[[1]][#], solFinal[[2]][#], solFinal[[4]][#]] & /@
    Subdivide[r0, rmax, 1000]], 16]];

Plot[
 Evaluate[solFinal[[1]][r]/v],
 {r, r0, 12},
 PlotRange -> All,
 AxesLabel -> {"r", "\[Phi]/v"},
 PlotLabel -> "Scalar profile"
]
'''
path='/mnt/data/GB_shooting_fixed_v3.wl'
open(path,'w',encoding='utf-8').write(code)
print(path, os.path.getsize(path))
print(code[:1000])

python
code = open('/mnt/data/GB_shooting_fixed_v3.wl', encoding='utf-8').read()
insert = r'''
rhsPhi[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, nn_?NumericQ, ll_?NumericQ] :=
  rhsVals[rr, pp, zz, nn, ll][[1]];
rhsF[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, nn_?NumericQ, ll_?NumericQ] :=
  rhsVals[rr, pp, zz, nn, ll][[2]];
rhsNu[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, nn_?NumericQ, ll_?NumericQ] :=
  rhsVals[rr, pp, zz, nn, ll][[3]];
rhsLam[rr_?NumericQ, pp_?NumericQ, zz_?NumericQ, nn_?NumericQ, ll_?NumericQ] :=
  rhsVals[rr, pp, zz, nn, ll][[4]];
'''
# Insert after rhsVals block ends before validSolutionQ
code = code.replace('validSolutionQ[sol_] :=', insert+'\nvalidSolutionQ[sol_] :=')
# replace NDSolve equations
code = code.replace('phi\'[r] == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[1]],\n      F\'[r]   == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[2]],\n      nu\'[r]  == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[3]],\n      lam\'[r] == rhsVals[r, phi[r], F[r], nu[r], lam[r]][[4]],',
'''phi'[r] == rhsPhi[r, phi[r], F[r], nu[r], lam[r]],
      F'[r]   == rhsF[r, phi[r], F[r], nu[r], lam[r]],
      nu'[r]  == rhsNu[r, phi[r], F[r], nu[r], lam[r]],
      lam'[r] == rhsLam[r, phi[r], F[r], nu[r], lam[r]],''')
open('/mnt/data/GB_shooting_fixed_v3.wl','w',encoding='utf-8').write(code)
print('wrote', len(code))
# show segment
idx=code.find('rhsPhi')
print(code[idx:idx+1000])