Я решаю систему дифф. уравнений (24), (25) и (27) ...

python
import sympy as sp
x=sp.symbols('x')
expr=x-1.6-0.8/(1-0.8*x)*(1+1.6*x)
sp.solve(sp.together(expr), x)


[0.625 - 1.61535599791501*I, 0.625 + 1.61535599791501*I]

python
import numpy as np, sympy as sp
def roots_for(ph,ell=1,rS=1,mu=0.1,h=1):
    v=mu/h
    x=sp.symbols('x')
    V=lambda y: h**2/4*(y**2-v**2)**2
    Vp=lambda y: h**2*y*(y**2-v**2)
    b=(1-rS**2*V(ph))/(rS-8*ell**2*ph*x)
    expr=x-16*ell**2*ph/rS**3 - 8*b*ell**2*ph/rS**2*(1+16*ell**2*ph*x/rS) - (Vp(ph)+16*ell**2*ph*V(ph)/rS**2)/b
    pol=sp.Poly(sp.together(expr).as_numer_denom()[0],x)
    roots=[complex(r) for r in np.roots([float(c) for c in pol.all_coeffs()])]
    out=[]
    for rr in roots:
        if abs(rr.imag)<1e-8:
            xval=rr.real
            bval=(1-rS**2*V(ph))/(rS-8*ell**2*ph*xval)
            out.append((xval,bval))
    return out
for ph in [0,0.001,0.005,0.01,0.02,0.05,0.08,0.1]:
    print(ph, roots_for(ph))

python
import sympy as sp
x=sp.symbols('x')
expr = x-1.6-0.8/(1-0.8*x)*(1+1.6*x)
print(sp.factor(sp.together(expr).as_numer_denom()[0]))
print(sp.expand(sp.together(expr).as_numer_denom()[0]))
# multiply to integer
poly=sp.Poly(sp.together(expr).as_numer_denom()[0],x)
print(poly)
print([c for c in poly.all_coeffs()])
print([sp.nsimplify(c) for c in poly.all_coeffs()])
print(sp.factor(sp.together(expr).as_numer_denom()[0]*(-sp.Rational(5,4)))) # maybe

python
import mpmath as mp

rS=1.0
ell=1.0
mu=0.1
h=1.0
v=mu/h
def V(phi): return h*h/4*(phi*phi-v*v)**2
def Vp(phi): return h*h*phi*(phi*phi-v*v)
def b1(phi,d):
    return (1-rS*rS*V(phi))/(rS-8*ell*ell*phi*d)
def horizon(phi,d):
    bb=b1(phi,d)
    return d - 16*ell**2*phi/rS**3 - (8*ell**2*phi*bb/rS**2)*(1+16*ell**2*phi*d/rS) - (Vp(phi)+16*ell**2*phi*V(phi)/rS**2)/bb

for ph in [0,0.001,0.01,0.05,0.1,0.2]:
    roots=[]
    for guess in [-10,-1,-.1,0,.1,1,10]:
        try:
            rt=mp.findroot(lambda d: horizon(ph,d), guess, tol=1e-12, maxsteps=100)
            if abs(mp.im(rt))<1e-8:
                r=float(mp.re(rt))
                if not any(abs(r-x)<1e-6 for x in roots):
                    roots.append(r)
        except Exception as e:
            pass
    print(ph, roots, [(d,b1(ph,d),horizon(ph,d)) for d in roots])

python
import numpy as np, math
def real_roots_for_phi(ph):
    roots=[]
    for guess in np.linspace(-5,5,101):
        try:
            rt=mp.findroot(lambda d: horizon(ph,d), guess, tol=1e-10, maxsteps=50)
            if abs(mp.im(rt))<1e-7:
                r=float(mp.re(rt))
                if abs(horizon(ph,r))<1e-6 and all(abs(r-x)>1e-5 for x in roots):
                    roots.append(r)
        except: pass
    return sorted(roots)

phis=np.linspace(0,0.12,25)
for ph in phis:
    roots=real_roots_for_phi(ph)
    print(f"{ph:.4f}", roots[:3])

python
import numpy as np, math
from scipy.integrate import solve_ivp
from scipy.optimize import brentq, root_scalar
rS=1.0
ell=1.0
mu=0.1
h=1.0
v=mu/h
eps=1e-4
rmax=50.0
def V(phi): return h*h/4*(phi*phi-v*v)**2
def Vp(phi): return h*h*phi*(phi*phi-v*v)
def b1_val(phi,d):
    return (1-rS*rS*V(phi))/(rS-8*ell*ell*phi*d)
def horizon(phi,d):
    bb=b1_val(phi,d)
    return d - 16*ell**2*phi/rS**3 - (8*ell**2*phi*bb/rS**2)*(1+16*ell**2*phi*d/rS) - (Vp(phi)+16*ell**2*phi*V(phi)/rS**2)/bb
def dphi_roots(phi):
    # find all roots in a range, including singular maybe
    xs=np.linspace(-50,50,20001)
    vals=[]
    roots=[]
    prevx=None; prevv=None
    for x in xs:
        denom = rS-8*ell*ell*phi*x
        if abs(denom)<1e-8:
            prevx=None; prevv=None
            continue
        try:
            y=horizon(phi,x)
            if not np.isfinite(y):
                prevx=None; prevv=None
                continue
        except Exception:
            prevx=None; prevv=None
            continue
        if prevv is not None:
            if y==0 or prevv*y<0:
                try:
                    rt=brentq(lambda d: horizon(phi,d), prevx, x, maxiter=100)
                    if abs(horizon(phi,rt))<1e-6 and all(abs(rt-r)>1e-5 for r in roots):
                        roots.append(rt)
                except Exception as e:
                    pass
        prevx=x; prevv=y
    return roots
for ph in [0,0.005,0.01,0.02,0.03,0.035,0.037,0.038,0.04]:
    print(ph,dphi_roots(ph))

python
def dphi_small(phi):
    roots=dphi_roots(phi)
    roots=[r for r in roots if b1_val(phi,r)>0]
    if not roots: 
        raise ValueError("no root")
    # choose smaller positive root continuous with 24 phi
    return min(roots, key=lambda x:abs(x-24*ell**2*phi/rS**3))

def rhs(r,y, include_nu=False):
    # y [p,z,lam] or [p,z,nu,lam]
    p,z=y[0],y[1]
    l=y[3] if len(y)==4 else y[2]
    E=math.exp(-l)
    eL=math.exp(l)
    A=(1/r)*(1 - (8*ell**2/r)*(1-3*E)*p*z)
    nuP=((eL-1)/r**2 + 0.5*z*z - eL*V(p))/A
    lamP=(-(eL-1)/r**2 - (16*ell**2/r**2)*(1-E)*p*z + 0.5*z*z + eL*V(p))/A
    # linear system
    a

python
from scipy.optimize import minimize_scalar
for ph in [0.04,0.05,0.1]:
    res=minimize_scalar(lambda d: horizon(ph,d)**2 if abs(1-8*ph*d)>1e-8 else 1e10, bounds=(-5,5), method='bounded')
    print(ph,res.x,horizon(ph,res.x), b1_val(ph,res.x),res.fun)

python
def vals(r,y):
    p,z,nu,l=y
    E=math.exp(-l); eL=math.exp(l)
    A=(1/r)*(1 - (8*ell**2/r)*(1-3*E)*p*z)
    nuP=((eL-1)/r**2 + 0.5*z*z - eL*V(p))/A
    lamP=(-(eL-1)/r**2 - (16*ell**2/r**2)*(1-E)*p*z + 0.5*z*z + eL*V(p))/A
    a11=1 + (16*ell**2*E/r)*p*z
    a12=(16*ell**2*E/r)*nuP*p
    a21=(8*ell**2*p/r**2)*(1-E)
    det=a11-a12*a21
    return A,nuP,lamP,det,E,eL,p,z,l
for ph in [0.001,0.01,0.03]:
    sol=integrate(ph,rmax=20,method='Radau')
    print("\nph",ph,"end",sol.t[-1],"message",sol.message)
    print(vals(sol.t[-1], sol.y[:,-1]))
    # last few
    for i in [-5,-1]:
        print(sol.t[i], sol.y[:,i])

python
for ph in [1e-8,1e-6,1e-5,1e-4,5e-4,1e-3,2e-3]:
    try:
        sol=integrate(ph,rmax=50,method='Radau')
        status=sol.status
        end=sol.t[-1]
        phiend=sol.y[0,-1]
        print("phiS",ph,"end",end,"status",status,"phiend",phiend,"res", phiend-v)
    except Exception as e: print('err',ph,e)

python
for ph in np.linspace(0.00055,0.00095,9):
    sol=integrate(ph,rmax=50,method='Radau')
    print(f"{ph:.6f} end {sol.t[-1]:.3f} status {sol.status} phiend {sol.y[0,-1]:.5f} res {sol.y[0,-1]-v:.5f} msg {sol.message[:30]}")

python
for ph in np.linspace(0.000605,0.000645,9):
    sol=integrate(ph,rmax=50,method='Radau')
    print(f"{ph:.6f} end {sol.t[-1]:.3f} status {sol.status} phiend {sol.y[0,-1]:.5f} res {sol.y[0,-1]-v:.5f}")

python
def residual_good(ph, rmax=50):
    sol=integrate(ph,rmax=rmax,method='Radau')
    if sol.status!=0:
        return np.nan, sol
    return sol.y[0,-1]-v, sol
lo,hi=0.000610,0.000615
for _ in range(20):
    mid=(lo+hi)/2
    res, sol=residual_good(mid)
    print(mid,res, sol.status, sol.t[-1])
    if np.isnan(res) or res>0:
        hi=mid
    else:
        lo=mid
print("root bracket",lo,hi)
res,sol=residual_good((lo+hi)/2)
print(res,sol.y[:,-1], sol.t[-1])
print("d", dphi_small((lo+hi)/2), "b", b1_val((lo+hi)/2, dphi_small((lo+hi)/2)))