我有一段Python代码，它从文本文件中读取两列数据，并执行曲线拟合过程。代码基于numpy amd数学库，第一次运行代码时，我得到了一组结果。在此运行期间，机器还运行另一个Python代码来进行不同的分析。当我重新运行代码时，我发现得到了截然不同的结果。
为了确认这些结果的正确性，我在MATLAB中重新运行了计算，并且能够重现第二次运行Python代码获得的结果。任何关于为什么会发生这种情况的指示将不胜感激。如果有人需要测试，我可以私下分享数据和代码。

import  pandas  as  pd
import  numpy   as  np
import  typing  as  tp
import  math    as  mt
import  matplotlib.pyplot as plt

def curveFitting(energyData: tp.Dict[float, float], order: int) -> None:
    pesDataFrame    =   pd.DataFrame(columns = ["Distance", "Energy"])
    for key in energyData:
        pesDataFrame.loc[len(pesDataFrame)] = [key, energyData[key]]    #   Create a DataFrame to drop all energy values above 0 cm-1
    pesDataFrame    =   pesDataFrame[pesDataFrame["Energy"] <= 0]        #   Drop all energies greater than 0 cm-1
#    pesDataFrame    =   pesDataFrame[pesDataFrame["Distance"] > 1.4]
    distances       =   pesDataFrame["Distance"].to_numpy()
    energies        =   pesDataFrame["Energy"].to_numpy()
    rEquilibrium    =   distances[energies.argmin()]
    
    distances       =   np.delete(distances, energies.argmin())
    energies        =   np.delete(energies, energies.argmin())
    xRvalues        =   np.zeros(distances.size)
    print(rEquilibrium)

    for index, value in enumerate(distances):
        xRvalues[index] =   (pow(value, 3) - pow(rEquilibrium, 3)) / (pow(value, 3) + pow(rEquilibrium, 3)) # Calculate x(r) values according to Eq 14

    xMatrix         =   np.zeros((distances.size, order+1))
    for row, value in enumerate(xRvalues):
        for ordercount in range(0, order+1):
            xMatrix[row][ordercount] = pow(value, ordercount)           #   Calculate xMatrix according to Eq 19
    print(xMatrix.shape)
    
#    term1           =   np.zeros(distances.size)
    term2           =   np.zeros(distances.size)
    rescaledEnergy  =   np.zeros(distances.size)
    minEnergy       =   energies[energies.argmin()]
    dissocEnergy    =   0 - minEnergy               #   Calculate dissociation energy as difference between minimum and zero energy.
    for index, value in enumerate(energies):
        rescaledEnergy[index]   =   (value - minEnergy) / dissocEnergy
   
    for index, value in enumerate(distances):
#        term1[index]            =   np.emath.log(1 - mt.sqrt(rescaledEnergy[index]))    #   Eq 11
        term2[index]            =   np.emath.log(1 - mt.sqrt(rescaledEnergy[index]))    #   Eq 12
        
#        term1[index]            =   term1[index] / (rEquilibrium - value)
        term2[index]            =   term2[index] / (rEquilibrium - value)

    # First calculate beta from Eq 11
    beta, _, _, _   =   np.linalg.lstsq(xMatrix, term2, rcond = None)
    print(beta.shape)

    calculateBetas  =   np.dot(xMatrix, beta)
    errorVector     =   term2 - calculateBetas
    print(np.linalg.norm(errorVector))
#    calculateBetas      =   np.zeros(distances.shape)
#    for row in range(xMatrix.shape[0]):
#        calculateBetas[row] =   np.dot(xMatrix[row], beta)

    calculateEnergies1  =   np.zeros(distances.shape)
#    errorVector         =   np.zeros(distances.shape)

    for index, value in enumerate(distances):
        distance    =   value - rEquilibrium
        betaterm    =   mt.exp(-1*calculateBetas[index]*distance)
        nextterm    =   pow(1-betaterm, 2) * dissocEnergy
        calculateEnergies1[index]   =   nextterm + minEnergy
#        errorVector[index]          =   energies[index] - calculateEnergies1[index]
    error   =   np.linalg.norm(errorVector)
    plt.plot(distances, energies, "-r", label = "UCCSD(T) Energy")
    plt.plot(distances, calculateEnergies1, "-b", label = "Fitted Energies")
    plt.title("Fitting Order : {} and Error : {:10.5f}".format(order, error))
    plt.xlabel("Interatomic Distance")
    plt.ylabel("Energy (cm-1)")
    plt.legend()
    plt.savefig("uccsdt_fitting_order_{}.png".format(order), dpi = 600)

with open("energies_neutral_631g.dat") as energyFile:
    uhfEnergy       =   {}
    uccsdEnergy     =   {}
    uccsdtEnergy    =   {}

    for line in energyFile:
        line    =   line.split()
        dist    =   float(line[0])  #   Distance
        uhf     =   float(line[1])  #   UHF Energy
        uccsd   =   float(line[2])  #   UCCSD Energy
        uccsdt  =   float(line[3])  #   UCCSD Energy with UCCSD(T) corrections
        uhfEnergy[dist]     =   uhf
        uccsdEnergy[dist]   =   uccsd
        uccsdtEnergy[dist]  =   uccsdt

for order in range(0,20):
    curveFitting(energyData = uccsdtEnergy, order = order)
    plt.clf()

here提供了用于拟合的原始数据