2

我在这段代码中遇到数学域错误。但是相同的功能F在另一个程序（我没有使用fsolve）中确实有效，没有任何域错误。

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

我收到此错误：

    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
4

0 回答 0