Skip to content

Conversation

ChrisRackauckas
Copy link
Member

This update is required due to SciML/DiffEqProblemLibrary.jl#138. Its inherently just a numerical instability. https://github.com/SciML/DiffEqProblemLibrary.jl/pull/138/files#diff-416284a99fdfe358ae2a33c649aa3a715568ff3cf4d0d0cf62e6c7a183aa9fc7R111 this change is effectively random, looking at:

zz = out[1] + 3.0 * x[1] - 2.0 * x[1] * x[2] + 2.0 * x[1]^3
zzz = out[2] + x[2] - x[2]^2 - 1.0

out[1] += x[1] * (3.0 - 2.0 * x[2] + 2.0 * x[1]^2)
out[2] += x[2] * (1.0 - x[2]) - 1.0

if abs(out[1] - zz) > 1e5
    @show out[1] - zz, out[2] - zzz, x[1], x[2]
    @show 2.0 * x[2], 2.0 * x[1]^2
    @show 3.0 * x[1], 2.0 * x[1] * x[2], 2.0 * x[1]^3
    error()
end

I get:

(out[1] - zz, out[2] - zzz, x[1], x[2]) = (-2048.0, 0.0, 1.032129073975807e6, -1.109530258395812e6)
(2.0 * x[2], 2.0 * x[1] ^ 2) = (-2.219060516791624e6, 2.130580850692314e12)
(3.0 * x[1], 2.0 * x[1] * x[2], 2.0 * x[1] ^ 3) = (3.096387221927421e6, -2.2903568762924146e12, 2.1990344404556452e18)
(out[1] - zz, out[2] - zzz, x[1], x[2]) = (-4.503599627370496e15, 0.0, 2.027697193465798e10, -3.714249421601087e10)
(2.0 * x[2], 2.0 * x[1] ^ 2) = (-7.428498843202174e10, 8.223111816778149e20)
(3.0 * x[1], 2.0 * x[1] * x[2], 2.0 * x[1] ^ 3) = (6.0830915803973946e10, -1.5062746256024978e21, 1.6673980752436494e31)

with different thresholds. Basically when x[1] is large the error is floating point noise, and that floating point noise is enough for Broyden to be stable vs unstable.

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

This update is required due to SciML/DiffEqProblemLibrary.jl#138. Its inherently just a numerical instability. https://github.com/SciML/DiffEqProblemLibrary.jl/pull/138/files#diff-416284a99fdfe358ae2a33c649aa3a715568ff3cf4d0d0cf62e6c7a183aa9fc7R111 this change is effectively random, looking at:

```julia
zz = out[1] + 3.0 * x[1] - 2.0 * x[1] * x[2] + 2.0 * x[1]^3
zzz = out[2] + x[2] - x[2]^2 - 1.0

out[1] += x[1] * (3.0 - 2.0 * x[2] + 2.0 * x[1]^2)
out[2] += x[2] * (1.0 - x[2]) - 1.0

if abs(out[1] - zz) > 1e5
    @show out[1] - zz, out[2] - zzz, x[1], x[2]
    @show 2.0 * x[2], 2.0 * x[1]^2
    @show 3.0 * x[1], 2.0 * x[1] * x[2], 2.0 * x[1]^3
    error()
end
```

I get:


```
(out[1] - zz, out[2] - zzz, x[1], x[2]) = (-2048.0, 0.0, 1.032129073975807e6, -1.109530258395812e6)
(2.0 * x[2], 2.0 * x[1] ^ 2) = (-2.219060516791624e6, 2.130580850692314e12)
(3.0 * x[1], 2.0 * x[1] * x[2], 2.0 * x[1] ^ 3) = (3.096387221927421e6, -2.2903568762924146e12, 2.1990344404556452e18)
```

```
(out[1] - zz, out[2] - zzz, x[1], x[2]) = (-4.503599627370496e15, 0.0, 2.027697193465798e10, -3.714249421601087e10)
(2.0 * x[2], 2.0 * x[1] ^ 2) = (-7.428498843202174e10, 8.223111816778149e20)
(3.0 * x[1], 2.0 * x[1] * x[2], 2.0 * x[1] ^ 3) = (6.0830915803973946e10, -1.5062746256024978e21, 1.6673980752436494e31)
```

with different thresholds. Basically when `x[1]` is large the error is floating point noise, and that floating point noise is enough for Broyden to be stable vs unstable.
@ChrisRackauckas ChrisRackauckas merged commit d175f58 into master Apr 23, 2025
22 of 30 checks passed
@ChrisRackauckas ChrisRackauckas deleted the ChrisRackauckas-patch-1 branch April 23, 2025 10:59
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant