[SciPy-User] brentq solver gives a strange bug

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] brentq solver gives a strange bug

Yuxiang Wang
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

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: brentq solver gives a strange bug

Yuxiang Wang
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 <[hidden email]> 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
[hidden email]
+1 (434) 284-0836

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: brentq solver gives a strange bug

David Mikolas
Good catch! Sooner or later everyone gets bitten by roundoff. Printing or plotting at the first sign of trouble is a good reflex. The OP of my first stackexchange answer here http://stackoverflow.com/a/29482967/3904031 had run an algorithm which did not accommodate roundoff error for years without trouble. 

Just a side note, the SciPy brentq has unexpected behavior if nan appears in one of the limits. 


On Sat, Nov 12, 2016 at 1:15 AM, Yuxiang Wang <[hidden email]> wrote:
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 <[hidden email]> 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
[hidden email]
<a href="tel:%2B1%20%28434%29%20284-0836" value="+14342840836" target="_blank">+1 (434) 284-0836

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user