Mathematica code vs Python

[MiloChocolatito]

Hello, i have a Code in mathematica that works very well, the code is:

Output:
Block[{$MaxPrecision = 300, $MaxExtraPrecision = 300},
  (*WorkingPrecision\[Rule]MachinePrecision,
  Method\[Rule]{"GaussKronrodRule","SymbolicProcessing"\[Rule]0}*)
  nf = 1 (*(Subscript[n, f]/n)^(2/3) \[Congruent] Subscript[E, f]/
  Subscript[E, F] = 1*);
  \[Delta]e3D = 1/1000;
  G3D = 1/
    10000;
  
  
  Quiet[
   Print[gapmulist3D2e = Table[
      {
        N[T],
        N[SetPrecision[\[Mu], 20]],
        N[SetPrecision[Abs[\[CapitalDelta]e], 20]]
        } /.
       FindRoot[{
         G3D*NIntegrate[
            
            Sqrt[\[Epsilon]k]*(1/
               Sqrt[(\[Epsilon]k - \[Mu])^2 + (\[CapitalDelta]e)^2]*
               Tanh[Sqrt[(\[Epsilon]k - \[Mu])^2 + \
(\[CapitalDelta]e)^2]/(2*T)]),
            {\[Epsilon]k, nf^(2/3), nf^(2/3) + \[Delta]e3D},
            WorkingPrecision -> 20,
            Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0}
            ] == (nf^(2/3) + (\[Delta]e3D/2) - \[Mu]),
         
         
         1 == ((3*(\[CapitalDelta]e)^2)/(8*G3D))(*bosones 2eCPs con K=
           0, es decir,  2Subscript[n, 0](T) = 3\!\(
\*SubsuperscriptBox[\(\[CapitalDelta]\), \(e\), \(2\)]/8\)G*)
           + (6/2^(3/2))*NIntegrate[
             
             Sqrt[\[Epsilon]K]*(Coth[(
                 2*nf^(2/3) + \[Delta]e3D + \[Epsilon]K - 2*\[Mu])/(
                 2*T)] - 1),(*\[Epsilon](K) = \[Epsilon](0) + 
             Subscript[\[Epsilon], K] = 2Subsuperscript[E, f, 3/2] + 
             Subscript[\[HBar]\[Omega], D] + Subscript[\[Epsilon], K]
              relacion de dispersion*) 
             {\[Epsilon]K, 0, Infinity}
             ]
           + (3/4)*NIntegrate[
             
             Sqrt[\[Epsilon]k]*(1 - ((\[Epsilon]k - \[Mu])/
                  Abs[\[Epsilon]k - \[Mu]] Tanh[
                   Abs[\[Epsilon]k - \[Mu]]/(2*T)])),
             {\[Epsilon]k, 0, nf^(2/3)},
             WorkingPrecision -> 20
             ] 
           + (3/4)*NIntegrate[
             
             Sqrt[\[Epsilon]k]*(1 - ((\[Epsilon]k - \[Mu])/
                  Sqrt[(\[Epsilon]k - \[Mu])^2 + (\[CapitalDelta]e)^2]
                   Tanh[Sqrt[(\[Epsilon]k - \[Mu])^2 + \
(\[CapitalDelta]e)^2]/(2*T)])),
             {\[Epsilon]k, nf^(2/3), nf^(2/3) + \[Delta]e3D},
             WorkingPrecision -> 20
             ]
           + (3/4)*NIntegrate[
             
             Sqrt[\[Epsilon]k]*(1 - ((\[Epsilon]k - \[Mu])/
                  Abs[\[Epsilon]k - \[Mu]] Tanh[
                   Abs[\[Epsilon]k - \[Mu]]/(2*T)])),
             {\[Epsilon]k, nf^(2/3) + \[Delta]e3D, Infinity},
             WorkingPrecision -> 20
              ]
         },
        {
         {\[Mu], 1},
         {\[CapitalDelta]e, 1/100000}
         },
        Method -> {"Newton", "StepControl" -> "TrustRegion"},
        WorkingPrecision -> 18
        ],
      {T, 1/10000000, 1/100000, 1/10000000}
      ]
    ]
   ]
  ] // Timing

I tried to "translate" in python but i can't reproduce the results, this is my best try:

import numpy as np
from scipy.integrate import quad
from scipy.optimize import least_squares
import csv

np.set_printoptions(precision=18)


def solve_system(T, initial_guess):
    def eq1(mu, Delta, T):
        def integrand(epsilon_k, mu, Delta, T):
            energy = np.sqrt((epsilon_k - mu)**2 + Delta**2 + 1e-20)
            return epsilon_k**0.5 / energy * np.tanh(energy / (2 * max(T, 1e-20)))
        integral, _ = quad(integrand, 1, 1 + 1/1000, args=(mu, Delta, T), epsabs=1e-12, epsrel=1e-12, limit=1000)
        return 1 + (1/2000) - mu - (1/10000) * integral

    def eq2(mu, Delta, T):
        term1 = (3 * Delta**2) / (8 * 1e-4)
        
        def integrand2(epsilon_k, mu, T):
            arg = (2 + 1/1000 + epsilon_k - 2*mu) / (2 * max(T, 1e-20))
            return epsilon_k**0.5 * ((1 / np.tanh(arg + 1e-20)) - 1)
        term2, _ = quad(integrand2, 0, np.inf, args=(mu, T), epsabs=1e-12, epsrel=1e-12, limit=1000)
        term2 *= 6 / (2**(3/2))
        
        def integrand3(epsilon_k, mu, T):
            diff = epsilon_k - mu
            return epsilon_k**0.5 * (1 - np.sign(diff) * np.tanh(np.abs(diff) / (2 * max(T, 1e-20))))
        term3, _ = quad(integrand3, 0, 1, args=(mu, T), epsabs=1e-12, epsrel=1e-12, limit=1000)
        term3 *= 3/4

        def integrand4(epsilon_k, mu, Delta, T):
            energy = np.sqrt((epsilon_k - mu)**2 + Delta**2 + 1e-20)
            return epsilon_k**0.5 * (1 - (epsilon_k - mu)/energy * np.tanh(energy / (2 * max(T, 1e-20))))
        term4, _ = quad(integrand4, 1, 1 + 1/1000, args=(mu, Delta, T), epsabs=1e-12, epsrel=1e-12, limit=1000)
        term4 *= 3/4

        def integrand5(epsilon_k, mu, T):
            diff = epsilon_k - mu
            return epsilon_k**0.5 * (1 - np.sign(diff) * np.tanh(np.abs(diff) / (2 * max(T, 1e-20))))
        term5, _ = quad(integrand5, 1 + 1/1000, np.inf, args=(mu, T), epsabs=1e-12, epsrel=1e-12, limit=1000)
        term5 *= 3/4

        return term1 + term2 + term3 + term4 + term5 - 1

    def system(vars, T):
        mu, Delta = vars
        return [eq1(mu, Delta, T), eq2(mu, Delta, T)]

    sol = least_squares(
        lambda vars: system(vars, T),
        x0=initial_guess,
        bounds=([0.99999, 0], [1.1, 1e-6]),  # Cotas más amplias
        ftol=1e-10,
        xtol=1e-10
    )
    return sol.x


temperatures = np.arange(1e-7, 1e-5 + 1e-7, 1e-7)
results = []

current_guess = [1, 1/1000000]

for T in temperatures:
    mu, Delta = solve_system(T, current_guess)
    results.append({"T": T, "μ": mu, "Δ": Delta})
    current_guess = [mu, Delta]  
    print(f"T = {T:.1e}: μ = {mu:.15f}, Δ = {Delta:.15f}")

Any suggestion?

[Arush]

Hi, I see the problem because most sites won't do it for you it will take way too long for manual edit but I managed to do it here.

import time
import mpmath as mp

# set working precision (digits)
mp.mp.dps = 25

# physical parameters
nf = mp.mpf(1)
delta_e3D = mp.mpf('1e-3')
G3D = mp.mpf('1e-4')

# bounds in epsilon-space
eps_f = nf**(mp.mpf(2)/mp.mpf(3))
eps_f_top = eps_f + delta_e3D

def compute_gap_mu(T):
    """
    For a given temperature T, solve the two coupled equations for μ and Δe.
    Returns (mu, delta_e).
    """
    def system(vars):
        mu, de = vars

        # ---- first equation: gap equation ----
        def integrand1(ek):
            x = mp.sqrt((ek - mu)**2 + de**2)
            return mp.sqrt(ek) * (1/x * mp.tanh(x/(2*T)))
        I1 = mp.quad(integrand1, [eps_f, eps_f_top])

        eq1 = G3D * I1 - (eps_f + delta_e3D/2 - mu)

        # ---- second equation: number/conservation equation ----
        # term A: bosonic modes integral with coth-1
        def integrandA(eK):
            arg = (2*eps_f + delta_e3D + eK - 2*mu)/(2*T)
            return mp.sqrt(eK) * (mp.coth(arg) - 1)
        A = (6/mp.power(2,1.5)) * mp.quad(integrandA, [0, mp.inf])

        # term B1: k from 0 to eps_f
        def integrandB1(ek):
            x = ek - mu
            return mp.sqrt(ek) * (1 - mp.sign(x) * mp.tanh(abs(x)/(2*T)))
        B1 = (3/mp.mpf(4)) * mp.quad(integrandB1, [0, eps_f])

        # term B2: k from eps_f to eps_f_top (pairing shell)
        def integrandB2(ek):
            x = ek - mu
            denom = mp.sqrt(x**2 + de**2)
            return mp.sqrt(ek) * (1 - (x/denom)*mp.tanh(denom/(2*T)))
        B2 = (3/mp.mpf(4)) * mp.quad(integrandB2, [eps_f, eps_f_top])

        # term B3: k above eps_f_top
        def integrandB3(ek):
            x = ek - mu
            return mp.sqrt(ek) * (1 - mp.sign(x) * mp.tanh(abs(x)/(2*T)))
        B3 = (3/mp.mpf(4)) * mp.quad(integrandB3, [eps_f_top, mp.inf])

        eq2 = 1 - ((3*de**2)/(8*G3D) + A + B1 + B2 + B3)

        return [eq1, eq2]

    # initial guesses: μ=1, Δe=1e-5
    mu0 = mp.mpf(1)
    de0 = mp.mpf('1e-5')
    mu, de = mp.findroot(system, (mu0, de0), method='hybrid')
    return mu, de

def main():
    # build the list of (T, mu, |de|)
    results = []
    t_start = time.time()

    # loop T from 1e-7 to 1e-5 in steps of 1e-7
    T = mp.mpf('1e-7')
    T_end = mp.mpf('1e-5')
    dT = mp.mpf('1e-7')
    while T <= T_end + mp.mpf('1e-20'):
        mu, de = compute_gap_mu(T)
        results.append( (float(T), float(mu), float(abs(de))) )
        T += dT

    t_end = time.time()
    print(f"Timing: {t_end - t_start:.3f} s")
    # results is analogous to Mathematica's gapmulist3D2e
    for T, mu, de in results:
        print(f"T={T:.1e}, mu={mu:.6e}, |Δe|={de:.6e}")

if __name__ == "__main__":
    main()