如何使用Python解决非线性方程

Question

我有以下代码：

#!/usr/bin/env python
from scipy.optimize import fsolve
import math

h = 6.634e-27
k = 1.38e-16

freq1 = 88633.9360e6
freq2 = 88631.8473e6
freq3 = 88630.4157e6

def J(freq,T):
     return (h*freq/k)/(math.exp(h*freq/(k*T))-1)

def equations(x,y,z,w,a,b,c,d):
     f1 = a*(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
     f2 = b*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
     f3 = c*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))
     f4 = d*(J((freq3+freq1)/2,(y+w)/2)-J((freq3+freq1)/2,2.73))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
     return (f1,f2,f3,f4)

在上面的代码中，我定义了一些方程。这组方程有4个非线性方程，里面包含4个已知的变量a到d，还有4个未知数x、y、z和w。我想要以某种方式定义a到d，然后把它们输入到fsolve这个函数里，这样就能找到x、y、z和w的唯一解。这可能吗？

Answer 1

我没有遇到你说的“x未定义”的错误，但这可能是命令行出现了问题。不过在脚本中，这个操作是很简单的，不需要做任何修改。

Scipy的 fsolve 函数需要一个可以调用的函数作为输入。所以如果你想求解 f1，你可以这样做：

def f1(x, y, z, a):
    return a*(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))

root = fsolve(f1, x0=1, args=(1, 1, 1))

这个操作是为了求解方程 f1(x, y, z, a) = 0，x是我们要找的值，而y、z和a是额外的参数（在这个例子中，它们的值都是1）。所以根据你想求解的变量，那个变量应该放在参数的最前面，比如如果你想求解a，就需要写成 f1(a, x, y, z)。x0 是你对结果的初始估计。

最后，这个函数会返回一个包含结果的数组。

想了解更多信息，可以查看这个链接：http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

更新： 查看了关于 如何用Python求解一对非线性方程？ 的回答。

如果你想同时求解x、y、z和w，你需要把它们作为一个整体传入。所以你最终会得到这样的代码：

h = 6.634e-27
k = 1.38e-16

freq1 = 88633.9360e6
freq2 = 88631.8473e6
freq3 = 88630.4157e6

def J(freq,T):
    return (h*freq/k)/(math.exp(h*freq/(k*T))-1)

def equations(p, a, b, c, d):
    x, y, z, w = p
    f1 = a*(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
    f2 = b*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
    f3 = c*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))
    f4 = d*(J((freq3+freq1)/2,(y+w)/2)-J((freq3+freq1)/2,2.73))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
    return (f1,f2,f3,f4)

x, y, z, w = fsolve(equations, x0=(1, 1, 1, 1), args=(1, 1, 1, 1), xtol=1e-13) # x0=(x, y, z, w) args=(a, b, c, d)
print x, y, z, w # prints [2.73, 2.73, <some value>, 2.73]
print equations((x, y, z, w), 1, 1, 1, 1) # prints a tuple with values of order 1e-13 to 1e-16 (so we can assume they are indeed the roots)