Scipy-optimize.root_scalar 不一致地给出 f(a) 和 f(b) 必须有不同的符号

Question

我正在编写一些代码以将一种模型（“宏自旋模型”）拟合到另一种模型（“微磁模型”）。来自微磁模型的数据是从一个名为 mumax 的开源包生成的，并放入 pandas 数据帧以供访问。列标有从 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_scalar 来计算特定磁场值的磁化值。函数“fun”使用 macrospin_angle 来计算磁滞回线。最终，我将在 scipy 最小二乘拟合例程中使用“fun”。

我遇到的问题是 root_scalar 告诉我两个端点 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 @。 2*b*np.sin(theta) 987654328 @ 987654328和@ 987654329也有值0。

所以，当您尝试在括号内找到root [0, np.pi]时，微小的舍入误差可以使差异与达到合适的支架之间。您的测试值（B数组）都有8个小数;导致错误的值可能在不太有效数字中不同。

在任何情况下，括号的选择很糟糕，因为您最终可能会在@ 487654332 @或theta=np.pi。您应该设置@ 987654334或类似的东西。