Scipyoptimize.root_标量不一致地给出f（a）和f（b）必须有不同的符号

Question

我正在编写一些代码，将一种模型（“宏旋模型”）与另一种模型（“微磁模型”）相匹配。来自微磁模型的数据由一个名为mumax的开源软件包生成，并放入熊猫数据帧中进行访问。这些列标有0到90的角度（模拟中的磁场角度）和“B”或“M”（用于磁场和磁化）

对于每个角度，我提取特定模拟的B和M列，提取与正B值对应的数据，并将其用作分析宏销模型的输入

这是有关守则的一部分：

def macrospin_angle(Ku, Kg, phi, Ms, B):
    #This calculates the angle of a macrospin in a field oriented phi degrees
    # This works for a single value of B
    theta0 = np.pi/2
    b = B*Ms/(2*Ku + 2*Kg*np.square(np.sin(phi)))
    f = lambda theta: np.sin(2*(theta-theta0)) + 2*b*np.sin(theta)
    return opt.root_scalar(f, bracket = [0, np.pi]).root

def fun(B_vals, Ku, Kg, phi, Ms):
    theta = np.zeros((len(B_vals),1))
    M = np.zeros((len(B_vals),1))
    for idx, current_B in enumerate(B):
        theta[idx] = macrospin_angle(Ku, Kg, phi, Ms, current_B)
        M[idx] = Ms*np.cos(theta[idx])
    return M


for angle in angles:    
    #First, extract the particular loop
    current_B = Full_Table[angle + "_B"]
    current_M = Full_Table[angle + "_M"]
    current_DT = pd.DataFrame()
    current_DT = pd.concat([current_DT, current_B, current_M], axis = 1)
    #Now extract the positive field range as two tables
    idx_max = np.argmax(current_DT[angle + "_B"].to_numpy())
    pos_DT_1 = current_DT[:idx_max]
    pos_DT_1 = pos_DT_1[pos_DT_1[angle + "_B"] > 0]
    pos_DT_2 = current_DT[idx_max:]
    pos_DT_2 = pos_DT_2[pos_DT_2[angle + "_B"] > 0]
    pos_DT_2 = pos_DT_2.iloc[::-1]
    #Average the two branches to get rid of hysteresis
    B_vals = 0.5*(pos_DT_1[angle + "_B"].to_numpy() + pos_DT_2[angle + "_B"].to_numpy())
    M_vals = 0.5*(pos_DT_1[angle + "_M"].to_numpy() + pos_DT_2[angle + "_M"].to_numpy())
    foo = fun(B_vals, 0.6E6, 0.06E6, 0, 1000e3)

函数“macrospin_angle”使用scipy.optimize.root_标量计算特定磁场值的磁化值。函数“fun”使用宏旋转角来计算磁滞回线。最后，我将在一个简单的最小二乘拟合例程中使用“fun”

我遇到的问题是根标量告诉我两个端点f（a）和f（b）没有不同的符号。但是，当我看f（a）和f（b）的值时，它们有不同的符号。更奇怪的是，如果我只是将macrospin_angle和fun复制到自己的脚本中，并直接使用pandas表中的B值，那么脚本就可以正常工作：

import numpy as np
from matplotlib import pyplot as plt
from scipy import optimize as opt
def macrospin_angle(x, B):
    #This calculates the angle of a macrospin in a field oriented phi degrees
    # This works for a single value of B
    theta0 = np.pi/2
    b = B*x[3]/(2*x[0] + 2*x[1]*np.square(np.sin(x[2])))
    f = lambda theta: np.sin(2*(theta-theta0)) + 2*b*np.sin(theta)
    return opt.root_scalar(f, bracket = [0, np.pi]).root
B = [0.00999848, 0.02999543, 0.04999238, 0.06998934, 0.08998629,\
       0.10998324, 0.12998018, 0.14997715, 0.16997411, 0.18997107,\
       0.20996803, 0.22996497, 0.24996191, 0.26995886, 0.28995581,\
       0.30995279, 0.32994974, 0.34994669, 0.36994366, 0.38994062,\
       0.40993755, 0.4299345 , 0.44993147, 0.46992839, 0.48992537,\
       0.50992234, 0.52991928, 0.54991621, 0.56991315, 0.58991012,\
       0.60990707, 0.62990403, 0.64990101, 0.66989795, 0.6898949 ,\
       0.70989188, 0.72988883, 0.7498858 , 0.7698827 , 0.78987965,\
       0.8098767 , 0.8298736 , 0.8498705 , 0.8698675 , 0.88986445,\
       0.9098614 , 0.9298584 , 0.9498553 , 0.96985224, 0.98984927,\
       1.00984623, 1.0298432 , 1.04984005, 1.06983695, 1.08983395,\
       1.10983085, 1.12982785, 1.14982485, 1.16982175, 1.1898187 ,\
       1.2098157 , 1.22981265, 1.2498096 , 1.2698066 , 1.28980355,\
       1.3098005 , 1.3297975 , 1.3497944 , 1.3697913 , 1.3897883 ,\
       1.40978525, 1.4297822 , 1.4497792 , 1.46977615, 1.4897731 ]

x = [0.5e6, 0.05e6, 0, 1000e3]

theta = np.zeros((len(B),1))
m = np.zeros((len(B),1))
for idx, current_B in enumerate(B):
    theta[idx] = macrospin_angle(x,current_B)
    m[idx] = x[3]*np.cos(theta[idx])
plt.plot(B,m)

上面的代码工作正常，但完全相同。我很迷路，所以任何提示都非常感谢

Answer 1

在该功能中：

f = lambda theta: np.sin(2*(theta-theta0)) + 2*b*np.sin(theta)

术语sin(2*(theta-theta0))在θ中是周期性的，周期为pi。此外，通过选择theta0=np.pi/2，术语在theta==0和theta=np.pi处的值为0。术语2*b*np.sin(theta)在theta==0和theta==np.pi处的值也为0

因此，当您尝试在括号[0, np.pi]内查找根时，微小的舍入错误可能会导致获得合适括号与否的差异。你的测试值（数组）都是8位小数；导致错误的值可能在不太重要的数字中有所不同

在任何情况下，括号的选择都是错误的，因为最终可能会得到平凡的解决方案theta=0或theta=np.pi。您应该设置bracket=[1e-9, np.pi-1e-9]或类似的设置