CODEWITHSUNDEEP
pythonscipysolver
Olmo
Olmo

Reputation: 415

fsolve from scipy.optimize fails

I am trying to compute the roots of a vectorvalued functions, but the solution I get from fsolve is entirely wrong. I do get the warning: 175: RuntimeWarning: The iteration is not making good progress, as measured by the improvement from the last ten iterations. warnings.warn(msg, RuntimeWarning)

Here comes the working example:


import numpy as np 
from scipy.optimize import fsolve


def AMOC(amoc_state,
         gamma= 1/0.34,
         theta = 1,
         mu = 7.5,
         sigma = 0.85):
    T = amoc_state[0]
    S = amoc_state[1]
    dT = -gamma * (T-theta) - T * (1+ mu*np.abs(T-S)) 
    dS = sigma-S*(1+mu*np.abs(T-S))
    return (dT, dS) 

test = fsolve(AMOC, (0.8,0.3), 2.3,xtol = 1e-7, factor = 0.1, maxfev = 1000)

this gives me test = array([0.57372733, 0.47660354]) which is obviously no root of the AMOC function since

AMOC((0.57372733, 0.47660354), gamma = 2.3) returns (-0.01121948095040759, 0.026224895476309684)

any help is very much appreciated!

Upvotes: 0

Views: 1612

Answers (1)

joni
joni

Reputation: 7157

As already mentioned in the comments, the abs function is not differentiable and so is your function AMOC. This is important since fsolves underlying algorithm uses the derivative information in order to solve F(x) = 0 and thus, assumes that the function F is differentiable.

However, as a workaround you can replace abs(x) by the smooth approximation sqrt(x**2 + 1.0e-4):

def AMOC(amoc_state, gamma= 1/0.34, theta = 1, mu = 7.5, sigma = 0.85):
    T, S = amoc_state
    smooth_abs = np.sqrt((T-S)**2 + 1.0e-4)
    dT = -gamma * (T-theta) - T * (1+ mu*smooth_abs) 
    dS = sigma-S*(1+mu*smooth_abs)
    return np.array((dT, dS))

Next, you can try another initial guess as Warren already mentioned. Otherwise, you can solve the problem as equivalent minimization problem by means of scipy.optimize.minimize. In my experience, the solvers are a bit more robust:

from scipy.optimize import minimize

res = minimize(lambda x: np.sum(AMOC(x, 2.3)**2), x0=(0.8,0.3))

yields the solution array([2.25836418e-01, 1.38576534e-06]).

Upvotes: 1

Related Questions