how to avoid runtime errors with scipy.optimize.fsolve

Question

Question

I am trying to numerically solve a non linear system of algebraic equations using scipy.optimize.fsolve.

I solve the system for several different values of its parameters (k1, k2, k3 below). For some values of the parameters fsolve finds the correct solution, while for others the following warning is occurs

RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)

and from the result in these cases it is clear that something has gone wrong since h(result) is not a vector of zeros. There is no doubt about the fact that the solution does exist and it has nothing qualitatively different from those cases where the correct solution is found.

What is usually recommended in these cases? Is it an issue of initial conditions?

Examples

In the following I will show the idea of how I solve the system of equations:

import numpy as np
from scipy.optimize import fsolve

# Parameters for the system of equations
k1, k2, k3 = 2., 4.5, 0.1

# Function for evaluating the residual of the system
def h(x):
    x1, x2, x3=x
    eqn1 = k1*x1 + k2*x2**2 - x3
    eqn2 = x1 - 2.*x2 + k3*x3 
    eqn3 = x1 + x2 + x3 - 1.
    return eqn1, eqn2, eqn3

# An initial guess
x0 = np.array([1., 0.5, 1.])

# Get the solution of the system of equations
result = fsolve(h, x0=x0, xtol=1e-5)

# Now, verify that the solution is correct
print(h(result)) # Should be near (0,0,0)

This works great sometimes, but for some values of k1, k2, k3 it raises the RuntimeWarning discussed above and returns the wrong result

# Bad parameters
k1, k2, k3 = 2., 4.5, -1.

# Bad result
result = fsolve(h, x0=x0, xtol=1e-5)

# Verification shows the residual is not near (0,0,0)
print(h(result))

You can (hopefully) solve this by tweaking the solving method or using a different one. — Stop harming Monica, Aug 15 '19 at 18:56
Numeric methods can often be parametrized. With `fsolve` you can choose the starting point, the way the jacobian is conputed and more. — Stop harming Monica, Aug 15 '19 at 19:01
TY @Nils I edited those two mistakes in my question. I was wondering more generally what I should do in order to find the solutions in cases like this. Is it a matter of changing the initial conditions? Are there better ways than `fsolve` ? — 3sm1r, Aug 15 '19 at 19:07
Typically will numeric methods give you bad/no results if the starting conditions are bad. So in general you will recieve better results by choosing better initial values. One way to achieve that automatically would be two call fsolve in a loop with varying initial parameters and checking which results solve your equation with a smaller derivation from 0. — Nils, Aug 15 '19 at 19:22
@3sm1r Could you share an example of `k1,k2,k3` where you run into the `RuntimeWarning` you mention? `k1,k2,k3=np.array([2.,4.5,0.1])` runs fine for me. — Moormanly, Aug 15 '19 at 20:33
@Moormanly Ty for your help. I have edited the question according to your suggestion. I just want to stress for clarity that the problem is not that the solution does not exist. I know for sure that it does but the algorithm does not find it. — 3sm1r, Aug 15 '19 at 21:14

Moormanly · Answer 1 · 2019-08-15T22:34:39.107

I don't know how to get rid of your RuntimeWarning using fsolve. If you don't mind solving your system of equations using some algebra instead, here's how you can do that.

Start with your original system of equations

def h(x):
    x1, x2, x3=x
    eqn1 = k1*x1 + k2*x2**2 - x3
    eqn2 = x1 - 2.*x2 + k3*x3 
    eqn3 = x1 + x2 + x3 - 1.
    return eqn1, eqn2, eqn3

matrixify it

mat = np.array([
    [ k1, 0, k2, -1  ],
    [ 1, -2, 0,   k3 ],
    [ 1,  1, 0,   1  ]
], dtype=float)

vec = np.array([0, 0, 1], dtype=float)

def h(x):
    x1, x2, x3=x
    return mat @ (x1, x2, x2**2, x3) - vec

find a 4d vector y solving mat @ y = vec as well as a vector ns in the nullspace of mat

y = np.linalg.lstsq(mat, vec, rcond=None)[0]

from scipy.linalg import null_space
ns = null_space(mat).flatten()

now use the fact that z = y + a * ns will also solve mat @ z = vec for any scalar a. This allows us to find some z satisfying z[1]**2 = z[2]. For such a z, we have x = z[[0,1,3]] as the solution to the original system of equations.

Plugging z = y + a * ns into the constraint z[1]**2 = z[2], we need (y[1] + a * ns[1])**2 = y[2] + a * ns[2]. Solving this quadratic equation for a yields

a = np.roots([ns[1]**2, 2 * y[1] * ns[1] - ns[2], y[1]**2 - y[2]])

Since you claim the original system of equations had a solution, a should hold real values.

# This assertion will fail if (and only if) the original system has no solutions.
assert np.isreal(a).all()

Finally, we can get both solutions to the original system of equations

soln1 = (y + a[0] * ns)[[0,1,3]]
soln2 = (y + a[1] * ns)[[0,1,3]]

how to avoid runtime errors with scipy.optimize.fsolve

Question

Examples

1 Answers1