paulnorthrop / itp

The Interpolate, Truncate, Project (ITP) Root-Finding Algorithm
https://paulnorthrop.github.io/itp/
GNU General Public License v3.0
9 stars 2 forks source link

Possible error in Wikipedia's ITP article #2

Closed waruyama closed 2 years ago

waruyama commented 2 years ago

I just learned about ITP yesterday and tried an implementation that is based on Wikipedia's pseudocode implementation of the algorithm.

It says:

if yitp > 0 then b = xitp and yb = yitp
else if yitp < 0 then a = xitp and ya = yitp
else a = xitp and y = xitp

In your implementation you multiply yITP with the sign of yb in the conditions:

    if (yITP * signyb > 0.0) {
      new_b = xITP ;
      yb = yITP ;
    } else if (yITP * signyb < 0.0) {
      new_a = xITP ;
      ya = yITP ;
    } else {

I do not understand the algorithm, but changing the conditions to multiply yitp with the sign of yb (or simply yb) made the implementation work. Also, I know from other algorithms (like Regula Falsi) that this multiplication is done in the conditions.

Since you seem to have deep understanding of the ITP method, my question would be if you can confirm that there is an error in the algorithm in the Wikipedia article.

paulnorthrop commented 2 years ago

The pseudocode in the Wikipedia article "... assumes the initial values of ya and yb ... satisfy ya < 0 < yb ...", which is consistent with the plots to the right of this text. This assumption simplifies the code for updating the bracketing interval. More generally, we need to be able to handle cases where ya > 0 > yb, which is what the multiplication by the sign does.

waruyama commented 2 years ago

Thank you. My ya and yb did not satisfy this condition, that's why the implementation failed. So we can leave the Wikipedia article as it is and only do the multiplication in the code.