I need to solve a non-linear set of three equations using scipy.

However, I do not have any clue on which algorithm is suitable for my problem from a mathematical point of view (stability, convergence behaviour), since scipy provides a huge variety of different algorithms in the scipy.optimize module such as:

Here is my set of equations (updated version due to comments):

set of equations

Which algorithm is probably the best one for my problem?

  • $\begingroup$ What is $t_{k-1}$? Is there some sort of initial condition for the $Q(0)$? $\endgroup$ – Bill Barth Jun 21 '15 at 14:10
  • $\begingroup$ $t_{k-1}$ just symbolizes the time dependence of the temperatures $\vartheta_{foo}$ whereas $k-1$ means that this references to the value of the last calculation / discretization step resp. initial condition if we talk about the first step. I am not sure, but $Q(0)$ should be zero, since there is no heat transfered at init. $\endgroup$ – albert Jun 21 '15 at 14:29
  • $\begingroup$ Yeah, I get that, but your last equation shows $t_k=1$, so you should either remove it or provide something about the size of the time steps. $\endgroup$ – Bill Barth Jun 21 '15 at 14:31
  • $\begingroup$ @BillBarth: Since this is part of a self-written simulation software the step size depends on the number of discretisation elements specified by the user. In the simulation algorithm I converted the rotation of a wheel into time steps depending on the total number of discretisation points. Due to simplicity I assumed that the step width is equal to one second which is given if the wheel turns with 1 rpm and is represented by 360 discretisation segments / points. If you wish having a more physical approach I could add the appropiate units to all variables. $\endgroup$ – albert Jun 21 '15 at 14:41
  • $\begingroup$ OK, replacing $t_k$ by $\Delta t$ makes this a bit clearer. $\endgroup$ – Bill Barth Jun 21 '15 at 15:00

Since your problem is small, you're probably best off trying fsolve or root. Both of these are interfaces to MINPACK and call HYBRD or HYBRJ. Since calculating a Jacobian matrix for your system shouldn't be hard (either do it by hand, or use your favorite computer algebra system, like SymPy, Sage, Maple, or Mathematica), you should supply a Jacobian matrix. The HYBR functions in MINPACK use "Powell's hybrid method", which uses Newton's method and checks if Newton steps will be descent steps by comparing against least-squares minimization (specifically, does the Newton step also decrease the sum-of-squared-residuals). If Newton steps are not descent steps, then it falls back to the gradient of the sum-of-squared-residuals.

I think you're probably better off ignoring all of the alternatives for now. All of these other alternatives, broadly speaking, try to find low-dimensional structure in high-dimensional problems. Your problem has dimension 3, so the comparative savings of finding lower-dimensional structure aren't that great compared to, say, PDE discretizations that might have tens of thousands of variables.


scipy.optimize.root is a unified interface for the nonlinear solvers that you list in your third point, which are implemented in Python. On the other hand, scipy.optimize.fsolve is a wrapper around some Fortran functions (see docs). The interface is the same for both and I haven't seen any benchmarks (quick Google search), so just implement the function of the system and give them a go. Focus on correct implementation of this function and on choosing a sensible initial condition and then just try both root with its different methods and fsolve.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.