# [SciPy-User] brentq solver gives a strange bug

3 messages
Open this post in threaded view
|
Report Content as Inappropriate

## [SciPy-User] brentq solver gives a strange bug

 Dear all,When I run the following code snippet, everything works and x0 is within a and b:---import numpy as npfrom scipy.optimize import newton, brentqE = 73000A = 244.439436713B = 520.091701874n = 0.258964804689strain_pt = 0.00722080901839  # 0.0613590727341def eqn(stress_pt, strain_pt, A, B, n, E):    return stress_pt - A - B * (strain_pt - stress_pt / E) ** nx0, r = brentq(eqn, args=(strain_pt, A, B, n, E), a=A, b=strain_pt * E,               full_output=True)---But when I run the following complete code (actually not that long, and self-contained):---import numpy as npfrom scipy.optimize import brentqE = 73000A = 244.439436713B = 520.091701874n = 0.258964804689def eqn(stress_pt, strain_pt, A, B, n, E):    return stress_pt - A - B * (strain_pt - stress_pt / E) ** nstrain = np.logspace(0, 4, 100, endpoint=True) / 1e4stress = np.empty_like(strain)for i in range(strain.shape[0]):    if strain[i] <= A / E:        stress[i] = E * strain[i]    else:        stress[i], r = brentq(eqn, args=(strain[i], A, B, n, E), a=A,                              b=strain[i] * E, full_output=True)---Some solutions are out of the bracket!If I do:for i in range(strain.shape[0]):    if strain[i] <= A / E:        stress[i] = E * strain[i]    else:        stress[i], r = brentq(eqn, args=(strain[i], A, B, n, E), a=A,                              b=strain[i] * E, full_output=True)        if stress[i] < A:            print(r)I can see two data points being wrong. Obviously, at i == 46 and i == 69, we have x0 = 0.0 which is out of [a, b].Could anyone please check whether this is repeatable on your machine?Shawn _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user
Open this post in threaded view
|
Report Content as Inappropriate

## Re: brentq solver gives a strange bug

 Dear all,Found the error - as always, it turned out to be my problem in defining the function.When strain_pt - stress_pt / E are too close to zero, sometimes we have a round-off error to cause it to be negative (-1e-18 or so). Taking the power of that causes nan, which is why the solution failed.Thanks,ShawnOn Fri, Nov 11, 2016 at 1:06 AM, Yuxiang Wang wrote:Dear all,When I run the following code snippet, everything works and x0 is within a and b:---import numpy as npfrom scipy.optimize import newton, brentqE = 73000A = 244.439436713B = 520.091701874n = 0.258964804689strain_pt = 0.00722080901839  # 0.0613590727341def eqn(stress_pt, strain_pt, A, B, n, E):    return stress_pt - A - B * (strain_pt - stress_pt / E) ** nx0, r = brentq(eqn, args=(strain_pt, A, B, n, E), a=A, b=strain_pt * E,               full_output=True)---But when I run the following complete code (actually not that long, and self-contained):---import numpy as npfrom scipy.optimize import brentqE = 73000A = 244.439436713B = 520.091701874n = 0.258964804689def eqn(stress_pt, strain_pt, A, B, n, E):    return stress_pt - A - B * (strain_pt - stress_pt / E) ** nstrain = np.logspace(0, 4, 100, endpoint=True) / 1e4stress = np.empty_like(strain)for i in range(strain.shape[0]):    if strain[i] <= A / E:        stress[i] = E * strain[i]    else:        stress[i], r = brentq(eqn, args=(strain[i], A, B, n, E), a=A,                              b=strain[i] * E, full_output=True)---Some solutions are out of the bracket!If I do:for i in range(strain.shape[0]):    if strain[i] <= A / E:        stress[i] = E * strain[i]    else:        stress[i], r = brentq(eqn, args=(strain[i], A, B, n, E), a=A,                              b=strain[i] * E, full_output=True)        if stress[i] < A:            print(r)I can see two data points being wrong. Obviously, at i == 46 and i == 69, we have x0 = 0.0 which is out of [a, b].Could anyone please check whether this is repeatable on your machine?Shawn -- Yuxiang "Shawn" Wang, PhDBiomechanical experiments & simulations[hidden email]+1 (434) 284-0836 _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user