diff --git a/src/solve.cc b/src/solve.cc index e27d033b..195725f5 100644 --- a/src/solve.cc +++ b/src/solve.cc @@ -810,6 +810,28 @@ list_p Root::multiple_equation_solver(list_r eqs, list_r names, list_r guesses) } +static array_g jacobi_escape_singular(size_t n, + algebraic_r step, + list_r guesses) +// ---------------------------------------------------------------------------- +// Perturb guesses to escape a degenerate (singular-Jacobian) starting point +// ---------------------------------------------------------------------------- +// Returns the perturbed guesses array, or null on failure. +{ + for (size_t col = 0; col < n; col++) + { + algebraic_g pw = step * decimal::make(col + 1); + if (!pw || !rt.push(+pw)) + return nullptr; + } + array_g perturb = array::from_stack(n, 0); + array_g v = array_p(+guesses); + v = v + perturb; + record(jsolve, "Singular Jacobian, perturbed guesses to %t", +v); + return v; +} + + bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) // ---------------------------------------------------------------------------- // Compute a Jacobian when there is cross-talk between variables @@ -826,9 +848,10 @@ bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) int impr = Settings.SolverImprecision(); algebraic_g eps = algebraic::epsilon(impr); algebraic_g oeps = decimal::make(101,-2); - uint max = Settings.SolverIterations(); - uint iter = 0; - int errs = 0; + uint max = Settings.SolverIterations(); + uint iter = 0; + int errs = 0; + int singular_errs = 0; bool back = false; // Go backwards array_g j, v, d; algebraic_g magnitude, last, forward; @@ -845,20 +868,19 @@ bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) // Set all variables to current value of guesses list::iterator gi = guesses->begin(); + algebraic_g shuffle; + if (errs) + { + // Shuffle slightly around current position (same for all variables) + shuffle = pow(oeps, errs) + eps; + if (!shuffle) + goto error; + } for (object_p varo : *vars) { object_p valo = *gi; - if (errs) - { - if (algebraic_g valg = valo->as_algebraic()) - { - // Shuffle slightly around current position - algebraic_g pw = pow(oeps, errs) + eps; - if (!pw) - goto error; - valo = pw; - } - } + if (errs && valo->as_algebraic()) + valo = shuffle; if (!directory::store_here(varo, valo)) goto error; ++gi; @@ -979,13 +1001,27 @@ bool Root::jacobi_solver(list_g &eqs, list_g &vars, list_g &guesses) j = array::from_stack(n, n, true); v = array::from_stack(n, 0); d = v / j; - record(jsolve, "Jacobian %t values %t delta %t", +j, +v, +d); - v = array_p(+guesses); - v = v - d; + if (!d) + { + if (++singular_errs >= 5) + goto error; + rt.clear_error(); + algebraic_g step = pow(oeps, singular_errs); + if (!step) + goto error; + v = jacobi_escape_singular(n, step, guesses); + last = nullptr; + back = false; + } + else + { + singular_errs = 0; + record(jsolve, "Jacobian %t values %t delta %t", +j, +v, +d); + v = array_p(+guesses); + v = v - d; + } if (!v) goto error; - - // This is the new guesses guesses = +v; } // while (iter < max) diff --git a/src/tests.cc b/src/tests.cc index 9362a44f..e9019887 100644 --- a/src/tests.cc +++ b/src/tests.cc @@ -8073,6 +8073,13 @@ void tests::solver_testing() "{ X Y } { 0 0 } ROOT", ENTER) .expect("{ X=0.5 Y=0.86602 54037 84 }"); + step("Jacobian solver, circle and line with singular initial Jacobian") + .test(CLEAR, "{ 'X^2+Y^2=1' 'X+Y=0' } { X Y } { 0 0 } ROOT", ENTER) + .expect("{ X=-0.70710 67811 87 Y=0.70710 67811 87 }"); + + step("Jacobian solver, ln equations with singular Jacobian at initial guess") + .test(CLEAR, "{ 'LN(X)+Y=0' 'LN(Y)+X=0' } { X Y } { 1 1 } ROOT", ENTER) + .expect("{ X=0.56714 32904 1 Y=0.56714 32904 1 }"); step("Solving when the variable is initialized with a constant") .test(CLEAR, ("m=Ⓒme " "'MSlv(ⒺRelativityMassEnergy;[E];[1 eV])' "