数学描述
已知多组x、y, 求下面非线性方程的三个常系数k、b、f
y
(
x
)
=
1
−
(
e
−
k
x
)
(
b
(
1
−
e
−
x
b
+
f
(
x
−
b
(
1
−
e
−
x
b
)
)
)
y(x) = 1 - (e^{-kx})(b(1-e^{-\frac{x}{b}} + f(x - b(1 - e^{-\frac{x}{b}})))
y(x)=1−(e−kx)(b(1−e−bx+f(x−b(1−e−bx)))
将方程转换为标准的二元函数
h
(
x
,
y
∣
k
,
b
,
f
)
=
y
−
(
(
e
−
k
x
)
(
b
(
1
−
e
−
x
b
+
f
(
x
−
b
(
1
−
e
−
x
b
)
)
)
)
h(x, y | k, b, f) = y - ((e^{-kx})(b(1-e^{-\frac{x}{b}} + f(x - b(1 - e^{-\frac{x}{b}}))))
h(x,y∣k,b,f)=y−((e−kx)(b(1−e−bx+f(x−b(1−e−bx))))
形成如下三元非线性方程组
{
h
(
x
0
,
y
0
∣
k
,
b
,
f
)
=
0
h
(
x
1
,
y
1
∣
k
,
b
,
f
)
=
0
h
(
x
2
,
y
2
∣
k
,
b
,
f
)
=
0
\begin{cases} h(x_0, y_0|k, b, f) = 0 \\ h(x_1, y_1|k, b, f) = 0 \\ h(x_2, y_2|k, b, f) = 0 \\ \end{cases}
⎩⎪⎨⎪⎧h(x0,y0∣k,b,f)=0h(x1,y1∣k,b,f)=0h(x2,y2∣k,b,f)=0
代码求解
利用python的fsolve求解,代码如下:
import numpy as np
from scipy.optimize import fsolve, root
def func(X):
def h(x, y, k, b, f):
return y - (1 - np.exp(-k * x)) * ((b * (1 - np.exp(- x / b)) + f * (x - b * (1 - np.exp(- x / b)))))
k, b, f = X.tolist()
x = [8.833, 4.949, 20.005]
y = [4.864, 2.811, 8.126]
return [h(x[i], y[i], k, b, f) for i in range(len(x))]
res = fsolve(func, x0=[1, 1, 1])
print(res) # [0.31822147 5.13982664 0.20741606]
要点说明
-
fslove求得的是局部最优解, 所以不一定能得到全局最优解,所以迭代初始值x0对结果的影响可能很大, 不合适的初值可能导致结果不收敛或偏差太大,从而爆出如下警告
RuntimeWarning: The iteration is not making good progress, as measured by the improvement from the last ten iterations.
-
fsolve一般在使用时仅需提供待求解的函数func和迭代初始值x0
-
当求解的是方程组时, func要求返回方程列表(list或者np.array都可以);当求解的是单个方程时,直接返回方程就行
-
更多信息参见官网文档