[SciPy-User] brentq solver gives a strange bug
Yuxiang Wang
yw5aj at virginia.edu
Fri Nov 11 12:15:58 EST 2016
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,
Shawn
On Fri, Nov 11, 2016 at 1:06 AM, Yuxiang Wang <yw5aj at virginia.edu> wrote:
> Dear all,
>
> When I run the following code snippet, everything works and x0 is within a
> and b:
>
> ---
> import numpy as np
> from scipy.optimize import newton, brentq
>
>
> E = 73000
> A = 244.439436713
> B = 520.091701874
> n = 0.258964804689
> strain_pt = 0.00722080901839 # 0.0613590727341
>
>
> def eqn(stress_pt, strain_pt, A, B, n, E):
> return stress_pt - A - B * (strain_pt - stress_pt / E) ** n
>
>
> x0, 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 np
> from scipy.optimize import brentq
>
>
> E = 73000
> A = 244.439436713
> B = 520.091701874
> n = 0.258964804689
>
>
> def eqn(stress_pt, strain_pt, A, B, n, E):
> return stress_pt - A - B * (strain_pt - stress_pt / E) ** n
>
>
> strain = np.logspace(0, 4, 100, endpoint=True) / 1e4
> stress = 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, PhD
Biomechanical experiments & simulations
yw5aj at virginia.edu
+1 (434) 284-0836
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20161111/234e4645/attachment.html>
More information about the SciPy-User
mailing list