8

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 non-linear set (two linear, and a non-linear) of equations with 3 unknown variables:

$$ Q(t_k)=m_Sc_S\left(\vartheta_S(t_k)-\vartheta_S(t_{k-1})\right)\\ Q(t_k)=\dot{m}_Gc_{pG}\big(\vartheta_{GE}-\vartheta_{GA}(t_{k})\big)\Delta t\\ Q(t_k)=A\alpha\frac{\big(\vartheta_S(t_k)-\vartheta_{GE}\big)-\big(\vartheta_S(t_k)-\vartheta_{GA}(t_k)\big)}{\ln{\left(\frac{\vartheta_S(t_k)-\vartheta_{GE}}{\vartheta_S(t_k)-\vartheta_GA(t_k)}\right)}}\Delta t $$

where unknown variables are:

  • $\vartheta_S(t_k)$
  • $\vartheta_{GA}(t_k)$
  • $Q(t_k)$

and

  • $m_S = 6.868\,\mathrm{kg}$
  • $\dot{m}_G = 0.007\,\mathrm{kg/s}$
  • $c_S = 500\,\mathrm{J/kg K}$
  • $c_{pG} = 1005\,\mathrm{J/kg K}$
  • $\vartheta_S(0) = 273\,\mathrm{K}$
  • $\vartheta_{GA}(0) = 293\,\mathrm{K}$
  • $\vartheta_{GE} = 353\,\mathrm{K}$
  • $A = 3.733\,\mathrm{m}^2$
  • $\alpha = 99\,\mathrm{W/m}^2\mathrm{K}$
  • $\Delta t = t_k-t_{k-1} = 1\,\mathrm{s}$

Which algorithm is probably the best one for my problem?

Anton Menshov
  • 8,672
  • 7
  • 38
  • 94
albert
  • 223
  • 2
  • 7
  • What is $t_{k-1}$? Is there some sort of initial condition for the $Q(0)$? – Bill Barth Jun 21 '15 at 14:10
  • $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. – albert Jun 21 '15 at 14:29
  • 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. – Bill Barth Jun 21 '15 at 14:31
  • @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. – albert Jun 21 '15 at 14:41
  • OK, replacing $t_k$ by $\Delta t$ makes this a bit clearer. – Bill Barth Jun 21 '15 at 15:00

2 Answers2

7

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.

Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
  • This is probably a dumb question, but I am curious roughly speaking, how are the residuals in the least squares calculated? I understand how when fitting parameters to data, you have residuals between the model and exp, and you tweak parameters to minimize the sum of residuals... but in a X equations X unknowns, with no data to use to compare the model with, what are the residuals? – Charlie Crown Jan 08 '22 at 01:02
2

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.

astrojuanlu
  • 1,187
  • 2
  • 12
  • 25