Python Forum

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?