用 python 求解超越方程组
Solving a system of transcendental equations with python
假设我有以下四个等式:
- cos(x)/x = a
- cos(y)/y = b
- a + b = 1
- c sinc(x) = d sinc(y)
对于未知变量 x, y, a
和 b
。请注意 cos(x)/x=a
有多个解决方案。变量 y
也是如此。我只对 x
和 y
值感兴趣,它们是第一个正根(如果重要的话)。
您可以安全地假设 a, b, c
和 d
是已知的实常数,都是正数。
在 Mathematica 中,解决这个问题的代码类似于:
FindRoot[{Cos[x]/x == 0.2 a + 0.1,
Cos[y]/y == 0.2 b + 0.1,
a + b == 1.0,
1.03*Sinc[x] == Sinc[y]*1.02},
{{x, .1}, {y, .1}, {a, .3}, {b, .1}}]
因此 returns
{x -> 1.31636, y -> 1.29664, a -> 0.456034, b -> 0.543966}
虽然这很容易,但我不知道如何在 python 中做类似的事情。因此,如果有人可以指导我(或简单地告诉我如何)解决这个问题,我将不胜感激。
您可以使用 root
:
import numpy as np
from scipy.optimize import root
def your_funcs(X):
x, y, a, b = X
f = [np.cos(x) / x - 0.2 * a - 0.1,
np.cos(y) / y - 0.2 * b - 0.1,
a + b - 1,
1.03 * np.sinc(x) - 1.02 * np.sinc(y)]
return f
sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1])
print(sol2.x)
将打印
[ 1.30301572 1.30987969 0.51530547 0.48469453]
您的函数必须以计算结果为 0 的方式定义,例如a + b - 1
而不是 a + b = 1
.
快速检查:
print(your_funcs(sol2.x))
给予
[-1.9356960478944529e-11, 1.8931356482454476e-11, 0.0, -4.1039033282785908e-11]
所以,解法应该可以了(请注意e-11
基本为0)
或者,您也可以使用 fsolve
:
from scipy.optimize import fsolve
sol3 = fsolve(your_funcs, [0.1, 0.1, 0.3, 0.1])
结果相同:
[ 1.30301572 1.30987969 0.51530547 0.48469453]
您可以使用 args
参数传递额外的参数:
def your_funcs(X, fac_a, fac_b):
x, y, a, b = X
f = [np.cos(x) / x - fac_a * a - 0.1,
np.cos(y) / y - fac_b * b - 0.1,
a + b - 1,
1.03 * np.sinc(x) - 1.02 * np.sinc(y)]
return f
sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1], args=(0.2, 0.2))
print(sol2.x)
这给你 "old" 输出:
[ 1.30301572 1.30987969 0.51530547 0.48469453]
如果你运行
sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1], args=(0.4, 0.2))
print(sol2.x)
然后您将收到:
[ 1.26670224 1.27158794 0.34096159 0.65903841]
假设我有以下四个等式:
- cos(x)/x = a
- cos(y)/y = b
- a + b = 1
- c sinc(x) = d sinc(y)
对于未知变量 x, y, a
和 b
。请注意 cos(x)/x=a
有多个解决方案。变量 y
也是如此。我只对 x
和 y
值感兴趣,它们是第一个正根(如果重要的话)。
您可以安全地假设 a, b, c
和 d
是已知的实常数,都是正数。
在 Mathematica 中,解决这个问题的代码类似于:
FindRoot[{Cos[x]/x == 0.2 a + 0.1,
Cos[y]/y == 0.2 b + 0.1,
a + b == 1.0,
1.03*Sinc[x] == Sinc[y]*1.02},
{{x, .1}, {y, .1}, {a, .3}, {b, .1}}]
因此 returns
{x -> 1.31636, y -> 1.29664, a -> 0.456034, b -> 0.543966}
虽然这很容易,但我不知道如何在 python 中做类似的事情。因此,如果有人可以指导我(或简单地告诉我如何)解决这个问题,我将不胜感激。
您可以使用 root
:
import numpy as np
from scipy.optimize import root
def your_funcs(X):
x, y, a, b = X
f = [np.cos(x) / x - 0.2 * a - 0.1,
np.cos(y) / y - 0.2 * b - 0.1,
a + b - 1,
1.03 * np.sinc(x) - 1.02 * np.sinc(y)]
return f
sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1])
print(sol2.x)
将打印
[ 1.30301572 1.30987969 0.51530547 0.48469453]
您的函数必须以计算结果为 0 的方式定义,例如a + b - 1
而不是 a + b = 1
.
快速检查:
print(your_funcs(sol2.x))
给予
[-1.9356960478944529e-11, 1.8931356482454476e-11, 0.0, -4.1039033282785908e-11]
所以,解法应该可以了(请注意e-11
基本为0)
或者,您也可以使用 fsolve
:
from scipy.optimize import fsolve
sol3 = fsolve(your_funcs, [0.1, 0.1, 0.3, 0.1])
结果相同:
[ 1.30301572 1.30987969 0.51530547 0.48469453]
您可以使用 args
参数传递额外的参数:
def your_funcs(X, fac_a, fac_b):
x, y, a, b = X
f = [np.cos(x) / x - fac_a * a - 0.1,
np.cos(y) / y - fac_b * b - 0.1,
a + b - 1,
1.03 * np.sinc(x) - 1.02 * np.sinc(y)]
return f
sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1], args=(0.2, 0.2))
print(sol2.x)
这给你 "old" 输出:
[ 1.30301572 1.30987969 0.51530547 0.48469453]
如果你运行
sol2 = root(your_funcs, [0.1, 0.1, 0.3, 0.1], args=(0.4, 0.2))
print(sol2.x)
然后您将收到:
[ 1.26670224 1.27158794 0.34096159 0.65903841]