CODEWITHSUNDEEP
scipysolver
Amadej Kristjan Kocbek
Amadej Kristjan Kocbek

Reputation: 191

Solver (fsolve in python.scipy) not working

I have been trying to solve the following system of equations using a scipy solver:

    from scipy.optimize import fsolve
    import math
    import numpy as np

    S0 = 1000
    u  = 1
    d  = 1

    delta_t = 1
    r = 0
    psi_r = (np.exp(-r*delta_t))**(-1)
    fi_r = np.exp(-r*delta_t)
    fi_2r = np.exp(-2*r*delta_t)

    def equations(p):
        p_u,p_m,p_d = p
        return (p_u+p_m+p_d - 1, fi_r*(p_u*(S0+u) + p_m*S0 + p_u*(S0-d)) -S0,fi_2r*
(p_u*p_u*(S0+2*u) + 2*p_m*p_u*(S0+u) + (2*p_u*p_d+p_m*p_m)*S0 + 2*p_m*p_d*(S0-d) + p_d*p_d*(S0-2*d)) - S0)

    p_u,p_m,p_d = fsolve(equations,(0.3,0.5,0.4))

    print(equations((p_u,p_m,p_d)))

The problem is that, despite the first equation stating that the sum of my unknowns should be 1, it never gives a result that would satisfy this. What I get is unexpected numbers on the order of 10 to the -12, or sometimes even negative numbers, which I know cannot be the correct solution.

I know I have to try several initial guesses, but what concerns me is that none of the guesses so far has given me probabilities which sum up to 1.

Upvotes: 0

Views: 2307

Answers (1)

Stelios
Stelios

Reputation: 5521

In general, when a system has multiple solutions, the solution returned by a numerical solver depends on the initial solution estimate as well as the algorithm implemented. However, in general, it is difficult to "control" the solution properties returned by the solver.

In your case, you require a solution where all its elements lie within the interval [0,1]. However, the fsolve algorithm converges to a solution that does not satisfy this property. Unfortunately, fsolve does not allow for imposing any constraints on the solution it returns (as is also the case for any other numerical equation solver, to the best of my knowledge).

A workaround for imposing constraints on the solution is to formulate the equation solving problem as a constrained optimization problem. One straightforward approach is the following.

Suppose you need to solve the two equations [f(x)==0, g(x)==0] with respect to the vector x. Note that any solution of the system is also a minimizer of the function f(x)**2+g(x)**2, which may be obtained by calling scipy.optimize.least_squares. Now, scipy.optimize.least_squares accepts bounds for the solution it returns, which essentially allow you to find a solution of the original system of equation with the desired bound properties. 

from scipy.optimize import least_squares
sol = least_squares(equations,(0.3,0.5,0.4), bounds = (0,1))
print('solution: ', sol.x, '\n', 'equations:', equations(sol.x))

  solution:  [ 0.295   0.4101  0.295 ] 
  equations: (0.0, 0.0, 0.0)

Note that the returned solution is indeed a probability distribution and a root for the system of equations.

Upvotes: 2

Related Questions