python - 用 Python 求解联立方程和积分

Question

我正在尝试解决这些联立方程：

(1)

(2)

我认为答案是：

a = 0.0031
b = 0.6587

我试过了：

from sympy import Eq, solve
from sympy.abc import a,b
from scipy.integrate import quad

def integrand(theta):
    return 105*theta**2*(1-theta)**4

def integrate(a,b):
    return quad(integrand, a, b)[0]

sol = solve([ Eq(a*(1-a)**5, b*(1-b)**5),
              Eq(integrate(a,b), 0.95) ])

print(sol)
print({ s:sol[s].evalf() for s in sol })

但我得到了：

TypeError: cannot determine truth value of Relational

我希望计算机自动进行集成，但是我尝试手动集成：

sol = solve([ Eq(a*(1-a)**5 - b*(1-b)**5),
              Eq(  b**3/3 -b**4 +6*b**5/5 -2*b**6/3 +b**7/7
                 -(a**3/3 -a**4 +6*a**5/5 -2*a**6/3 +a**7/7)
                 -0.95/105) ])

print(sol)

但我得到了：

[]

所以我尝试了 nsolve，因为您可以设置求解器的起点：

from sympy import Symbol, nsolve

a = Symbol('a')
b = Symbol('b')

f1 = a*(1-a)**5 - b*(1-b)**5
f2 = b**3/3 -b**4 +6*b**5/5 -2*b**6/3 +b**7/7 - (a**3/3 -a**4 +6*a**5/5 -2*a**6/3 +a**7/7) -0.95/105

Sol = nsolve((f1,f2), (a,b), (0.5,0.5))
print(Sol)

但我得到了：

ZeroDivisionError: matrix is numerically singular

接下来我尝试为求解器提供我认为正确的答案（a0，b0）作为起点。

from sympy import Symbol, nsolve
from scipy.integrate import quad

a0 = 0.0031
b0 = 0.6587

a = Symbol('a')
b = Symbol('b')

f1 = a*(1-a)**5 - b*(1-b)**5
f2 = b**3/3 -b**4 +6*b**5/5 -2*b**6/3 +b**7/7 - (a**3/3 -a**4 +6*a**5/5 -2*a**6/3 +a**7/7) -0.95/105

Sol = nsolve((f1,f2), (a,b), (a0,b0))

a1 = Sol[0]
b1 = Sol[1]

print('a: ',a0,a1)
print('b: ',b0,b1)

def integrand(theta):
    return 105*theta**2*(1-theta)**4

def integrate(a,b):
    return quad(integrand, a, b)[0]

a,b=a0,b0
print(' ')
print('a,b=a0,b0')
print('a: ',a)
print('b: ',b)
print('a*(1-a)**5: ',a*(1-a)**5)
print('b*(1-b)**5: ',b*(1-b)**5)
I = integrate(a,b)
print('I: ',I)

a,b=a1,b1
print(' ')
print('a,b=a1,b1')
print('a: ',a)
print('b: ',b)
print('a*(1-a)**5: ',a*(1-a)**5)
print('b*(1-b)**5: ',b*(1-b)**5)
I = integrate(a,b)
print('I: ',I)

我有：

a:  0.00309652149716206
b:  0.658740229906034
a*(1-a)**5:  0.00304887526056235
b*(1-b)**5:  0.00304887526056235
I:  0.95

(a,b) 满足方程 1（它们有a*(1-a)**5 = b*(1-b)**5)方程 2（I 为 0.95）。正如我所说，我希望计算机自动进行积分：）

标签： pythonmathequationsimultaneous