fix user's guide to use new solver functions, fixes #712

Author: avehtari

src/bibtex/all.bib

@@ -1834,3 +1834,14 @@ @article{Vehtari+etal:2021:Rhat
   volume=16,
  pages={667--718}
 }
+
+@article{Timonen+etal:2023:ODE-PSIS,
+  title={An importance sampling approach for reliable and efficient inference in {Bayesian} ordinary differential equation models},
+  author={Timonen, Juho and Siccha, Nikolas and Bales, Ben and L{\"a}hdesm{\"a}ki, Harri and Vehtari, Aki},
+  journal={Stat},
+  year={2023},
+  volume = 12,
+  number = 1,
+  pages = {e614} 
+}
+

src/stan-users-guide/algebraic-equations.qmd

@@ -84,7 +84,7 @@ vector for parameters if the system does not involve data or parameters.
 Let's suppose $\theta = (3, 6)$. To call the algebraic solver, we need to
 provide an initial guess. This varies on a case-by-case basis, but in general
 a good guess will speed up the solver and, in pathological cases, even determine
-whether the solver converges or not. If the solver does not converge, the metropolis
+whether the solver converges or not. If the solver does not converge, the Metropolis
 proposal gets rejected and a warning message, stating no acceptable solution was
 found, is issued.
 
@@ -107,7 +107,7 @@ transformed parameters {
   vector[2] theta = [3, 6]';
   vector[2] y;
 
-  y = algebra_solver_newton(system, y_guess, theta, x_r, x_i);
+  y = solve_newton(system, y_guess, theta, x_r, x_i);
 }
 ```
 
@@ -137,24 +137,26 @@ For instance, it might make "physical sense" for a solution to be positive or ne
 
 On the other hand, a system may not have a solution (for a given point in the parameter
 space). In that case, the solver will not converge to a solution. When the solver fails to
-do so, the current metropolis proposal gets rejected.
+do so, the current Metropolis proposal gets rejected.
 
 ## Control parameters for the algebraic solver {#algebra-control.section}
 
-The call to the algebraic solver shown previously uses the default control settings. The solver
-allows three additional parameters, all of which must be supplied if any of them is
-supplied.
+The call to the algebraic solver shown previously uses the default control settings. The `_tol` variant of the solver function
+allows three additional parameters, all of which must be supplied.
 
 ```stan
-y = algebra_solver_newton(system, y_guess, theta, x_r, x_i,
-                          rel_tol, f_tol, max_steps);
+y = solve_newton_tol(system, y_guess, theta, x_r, x_i,
+                     scaling_step, f_tol, max_steps);
 ```
 
-The three control arguments are relative tolerance, function tolerance, and maximum
-number of steps. Both tolerances need to be satisfied. If one of them is not met, the
-metropolis proposal gets rejected with a warning message explaining which criterion
-was not satisfied. The default values for the control arguments are respectively
-`rel_tol = 1e-10` ($10^{-10}$), `f_tol = 1e-6` ($10^{-6}$), and `max_steps = 1e3` ($10^3$).
+For the Newton solver the three control arguments are scaling step, function tolerance, and maximum number of steps. For the Powell's hybrid method the three control arguments are relative tolerance, function tolerance, and maximum number of steps. If a Newton step is smaller than the scaling step tolerance, the code breaks, assuming the solver is no longer making significant progress. If set to 0, this constraint is ignored. For Powell's hybrid method the relative tolerance is the estimated relative error of the solver and serves to test if a satisfactory solution has been found. After convergence of the either solver, the proposed solution
+is plugged into the algebraic system and its norm is compared to the function tolerance. If the norm is below the function tolerance, the solution is deemed acceptable.  If the solver solver reaches the maximum number of steps, it stops and returns an error message. If one of the criteria is not met, the
+Metropolis proposal gets rejected with a warning message explaining which criterion
+was not satisfied. 
+
+
+The default values for the control arguments are respectively
+`scaling_step = 1e-3` ($10^{-3}$), `rel_tol = 1e-10` ($10^{-10}$), `f_tol = 1e-6` ($10^{-6}$), and `max_steps = 200` ($200$).
 
 ### Tolerance {-}
 
@@ -172,12 +174,12 @@ Smaller relative tolerances produce more accurate solutions but require more com
 #### Sensitivity analysis {-}
 
 The tolerances should be set low enough that setting them lower does not change the
-statistical properties of posterior samples generated by the Stan program.
+statistical properties of posterior samples generated by the Stan program. The sensitivity can be analysed using importance sampling without need to re-run MCMC with different tolerances as shown by @Timonen+etal:2023:ODE-PSIS.
 
 ### Maximum number of steps {-}
 
 The maximum number of steps can be used to stop a runaway simulation. This can arise in
 MCMC when a bad jump is taken, particularly during warmup. If the limit is hit, the
-current metropolis proposal gets rejected. Users will see a  warning message stating the
+current Metropolis proposal gets rejected. Users will see a  warning message stating the
 maximum number of steps has been exceeded.