2

Im getting math domain error in this code. But the same function F does work in another program (I didn't use fsolve) without any domain error.

import numpy as np
from scipy.integrate import quad
from scipy.optimize import fsolve
from math import sqrt

def F(x, g, a, d, A, M):
  return A*(x**(1-g))/(1-g)+M*((x**(a+1))/(a+1)+(x**(d+1))/(d+1))

def C(q):
  return q+1

def intf1(x, rho, g, a, d, A, M):
  return (1/sqrt(F(rho, g, a, d, A, M)-F(x, g, a, d, A, M)))

def eq_to_solve(q, rho, g, a, d, A, M):
  return (quad(intf1,0,rho,args=(rho,g, a, d, A, M))[0]-quad(intf1,q,rho,args=(rho,g, a, d, A, M))[0]-(C(q)*q)/intf1(q, rho, g, a, d, A, M))

A = -4
M = 20
g = 0.6
a = 0.7
d = 2.1
rhomin = 10**(-6)
rhomax = 50
num_rho_values = 60
rholist = np.linspace(rhomin, rhomax, num=num_rho_values, endpoint=False)
lambdalist=[]
found_rholist = []
for rho in rholist:
    q_guess_list = np.linspace(1, rho, num = 10, endpoint=False)
    for q_guess in q_guess_list:
        q_val = fsolve(eq_to_solve, q_guess, args=(rho, g, a, d, A, M))
        if q_val > 0 and abs(eq_to_solve(q_val, rho, g, a, d, A, M)) < 10**(-6):
            lam = (C(q_val)*q)**2/(2*intf1(q_val, rho, g, a, d, A, M))
            print("found q= ", q_val)
            if lam > 0:
                print('found lam = ', lam)
                lambdalist.append(lam)
                found_rholist.append(rho)
            else:
                None
        else:
            None

I'm getting this error:

    62                                 q_guess_list = np.linspace(1, rho, num = 10, endpoint=False)
     63                                 for q_guess in q_guess_list:
---> 64                                         q_val = fsolve(eq_to_solve, q_guess, args=(rho, g, a, d, A, M))
     65                                         if q_val > 0 and abs(eq_to_solve(q_val, rho, g, a, d, A, M)) < 10**(-6):
     66                                                 lam = (C(q_val)*q)**2/(2*intf1(q_val, rho, g, a, d, A, M))

      6 
      7 def intf1(x, rho, g, a, d, A, M):
----> 8   return (1/sqrt(F(rho, g, a, d, A, M)-F(x, g, a, d, A, M)))
      9 
     10 def eq_to_solve(q, rho, g, a, d, A, M):

ValueError: math domain error
JohanC
  • 71,591
  • 8
  • 33
  • 66
mathfun
  • 101
  • 7
  • It would be helpful to have values for `rhomin`, `rhomax` and `num_rho_values`. Also, the indentation is off. – yad Jul 31 '20 at 01:02
  • Than you for pointing out those massing values. I cant use tab to indent in this editor, so I used spaces. – mathfun Jul 31 '20 at 01:10
  • You can use ``` before and after your code to paste it all in nicely. – JimmyCarlos Jul 31 '20 at 01:13
  • 2
    The domain error might indicate that you are trying to square root negative numbers. Negative numbers are quite likely here. Also you might want to avoid dividing by zero. Maybe `return (1/sqrt(1e-30 - abs(F(rho, ...)-F(x, ...))))`?? – JohanC Jul 31 '20 at 09:26
  • The example is too complicated. Boil it down to only a few lines of code such that the error still reproduces. This is a bit of work, but often makes the issue very visible. – Nico Schlömer Aug 03 '20 at 21:16

0 Answers0