如何使用Python求解一对非线性方程?

发布于 2024-12-25 03:36:05 字数 181 浏览 1 评论 0原文

使用 Python 求解一对非线性方程的(最佳)方法是什么。(Numpy、Scipy 或 Sympy)

例如:

  • x+y^2 = 4
  • e^x+ xy = 3

解决上述问题的代码片段会很棒

What's the (best) way to solve a pair of non linear equations using Python. (Numpy, Scipy or Sympy)

eg:

  • x+y^2 = 4
  • e^x+ xy = 3

A code snippet which solves the above pair will be great

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(9

起风了 2025-01-01 03:36:05

对于数值求解,您可以使用 fsolve:

http://docs.scipy.org/doc/scipy/reference/ generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve

from scipy.optimize import fsolve
import math

def equations(p):
    x, y = p
    return (x+y**2-4, math.exp(x) + x*y - 3)

x, y =  fsolve(equations, (1, 1))

print equations((x, y))

for numerical solution, you can use fsolve:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve

from scipy.optimize import fsolve
import math

def equations(p):
    x, y = p
    return (x+y**2-4, math.exp(x) + x*y - 3)

x, y =  fsolve(equations, (1, 1))

print equations((x, y))
硪扪都還晓 2025-01-01 03:36:05

简短回答:使用 fsolve

正如其他答案中提到的,针对您提出的特定问题,最简单的解决方案是使用类似 fsolve 的东西:

from scipy.optimize import fsolve
from math import exp

def equations(vars):
    x, y = vars
    eq1 = x+y**2-4
    eq2 = exp(x) + x*y - 3
    return [eq1, eq2]

x, y =  fsolve(equations, (1, 1))

print(x, y)

输出:

0.6203445234801195 1.8383839306750887

解析解?

你说如何“解决”,但有不同种类的解决方案。既然你提到了 SymPy,我应该指出这可能意味着分析数字解决方案之间的最大区别。您给出的特定示例没有(简单的)解析解,但其他非线性方程组有。当有现成的分析解时,SymPY 通常可以为您找到它们:

from sympy import *

x, y = symbols('x, y')
eq1 = Eq(x+y**2, 4)
eq2 = Eq(x**2 + y, 4)

sol = solve([eq1, eq2], [x, y])

输出:

⎡⎛ ⎛  5   √17⎞ ⎛3   √17⎞    √17   1⎞  ⎛ ⎛  5   √17⎞ ⎛3   √17⎞    1   √17⎞  ⎛ ⎛  3   √13⎞ ⎛√13   5⎞  1   √13⎞  ⎛ ⎛5   √13⎞ ⎛  √13   3⎞  1   √13⎞⎤
⎢⎜-⎜- ─ - ───⎟⋅⎜─ - ───⎟, - ─── - ─⎟, ⎜-⎜- ─ + ───⎟⋅⎜─ + ───⎟, - ─ + ───⎟, ⎜-⎜- ─ + ───⎟⋅⎜─── + ─⎟, ─ + ───⎟, ⎜-⎜─ - ───⎟⋅⎜- ─── - ─⎟, ─ - ───⎟⎥
⎣⎝ ⎝  2    2 ⎠ ⎝2    2 ⎠     2    2⎠  ⎝ ⎝  2    2 ⎠ ⎝2    2 ⎠    2    2 ⎠  ⎝ ⎝  2    2 ⎠ ⎝ 2    2⎠  2    2 ⎠  ⎝ ⎝2    2 ⎠ ⎝   2    2⎠  2    2 ⎠⎦

请注意,在此示例中,SymPy 找到所有解,并且不需要给出初始估计。

您可以使用 evalf 对这些解进行数值评估:

soln = [tuple(v.evalf() for v in s) for s in sol]
[(-2.56155281280883, -2.56155281280883), (1.56155281280883, 1.56155281280883), (-1.30277563773199, 2.30277563773199), (2.30277563773199, -1.30277563773199)]

数值解的精度

但是,大多数非线性方程组都没有合适的解析解,因此使用上面的 SymPy 非常有用,但并不普遍适用。这就是为什么我们最终寻找数值解,即使有数值解:
1)我们不能保证我们已经找到了所有解决方案,或者当有很多解决方案时,我们不能保证找到“正确”的解决方案。
2)我们必须提供一个初步的猜测,这并不总是那么容易。

接受我们想要数字解决方案后,像 fsolve 这样的东西通常可以满足您的所有需求。对于此类问题,SymPy 可能会慢得多,但它可以提供其他功能,更精确地找到(数字)解决方案:

from sympy import *

x, y = symbols('x, y')
nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1])
⎡0.620344523485226⎤
⎢                 ⎥
⎣1.83838393066159 ⎦

具有更高的精度:

nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1], prec=50)
⎡0.62034452348522585617392716579154399314071550594401⎤
⎢                                                    ⎥
⎣ 1.838383930661594459049793153371142549403114879699 ⎦

Short answer: use fsolve

As mentioned in other answers the simplest solution to the particular problem you have posed is to use something like fsolve:

from scipy.optimize import fsolve
from math import exp

def equations(vars):
    x, y = vars
    eq1 = x+y**2-4
    eq2 = exp(x) + x*y - 3
    return [eq1, eq2]

x, y =  fsolve(equations, (1, 1))

print(x, y)

Output:

0.6203445234801195 1.8383839306750887

Analytic solutions?

You say how to "solve" but there are different kinds of solution. Since you mention SymPy I should point out the biggest difference between what this could mean which is between analytic and numeric solutions. The particular example you have given is one that does not have an (easy) analytic solution but other systems of nonlinear equations do. When there are readily available analytic solutions SymPY can often find them for you:

from sympy import *

x, y = symbols('x, y')
eq1 = Eq(x+y**2, 4)
eq2 = Eq(x**2 + y, 4)

sol = solve([eq1, eq2], [x, y])

Output:

⎡⎛ ⎛  5   √17⎞ ⎛3   √17⎞    √17   1⎞  ⎛ ⎛  5   √17⎞ ⎛3   √17⎞    1   √17⎞  ⎛ ⎛  3   √13⎞ ⎛√13   5⎞  1   √13⎞  ⎛ ⎛5   √13⎞ ⎛  √13   3⎞  1   √13⎞⎤
⎢⎜-⎜- ─ - ───⎟⋅⎜─ - ───⎟, - ─── - ─⎟, ⎜-⎜- ─ + ───⎟⋅⎜─ + ───⎟, - ─ + ───⎟, ⎜-⎜- ─ + ───⎟⋅⎜─── + ─⎟, ─ + ───⎟, ⎜-⎜─ - ───⎟⋅⎜- ─── - ─⎟, ─ - ───⎟⎥
⎣⎝ ⎝  2    2 ⎠ ⎝2    2 ⎠     2    2⎠  ⎝ ⎝  2    2 ⎠ ⎝2    2 ⎠    2    2 ⎠  ⎝ ⎝  2    2 ⎠ ⎝ 2    2⎠  2    2 ⎠  ⎝ ⎝2    2 ⎠ ⎝   2    2⎠  2    2 ⎠⎦

Note that in this example SymPy finds all solutions and does not need to be given an initial estimate.

You can evaluate these solutions numerically with evalf:

soln = [tuple(v.evalf() for v in s) for s in sol]
[(-2.56155281280883, -2.56155281280883), (1.56155281280883, 1.56155281280883), (-1.30277563773199, 2.30277563773199), (2.30277563773199, -1.30277563773199)]

Precision of numeric solutions

However most systems of nonlinear equations will not have a suitable analytic solution so using SymPy as above is great when it works but not generally applicable. That is why we end up looking for numeric solutions even though with numeric solutions:
1) We have no guarantee that we have found all solutions or the "right" solution when there are many.
2) We have to provide an initial guess which isn't always easy.

Having accepted that we want numeric solutions something like fsolve will normally do all you need. For this kind of problem SymPy will probably be much slower but it can offer something else which is finding the (numeric) solutions more precisely:

from sympy import *

x, y = symbols('x, y')
nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1])
⎡0.620344523485226⎤
⎢                 ⎥
⎣1.83838393066159 ⎦

With greater precision:

nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1], prec=50)
⎡0.62034452348522585617392716579154399314071550594401⎤
⎢                                                    ⎥
⎣ 1.838383930661594459049793153371142549403114879699 ⎦
清欢 2025-01-01 03:36:05

如果您更喜欢 sympy,可以使用 nsolve

>>> nsolve([x+y**2-4, exp(x)+x*y-3], [x, y], [1, 1])
[0.620344523485226]
[1.83838393066159]

第一个参数是方程列表,第二个参数是变量列表,第三个参数是初始猜测。

If you prefer sympy you can use nsolve.

>>> nsolve([x+y**2-4, exp(x)+x*y-3], [x, y], [1, 1])
[0.620344523485226]
[1.83838393066159]

The first argument is a list of equations, the second is list of variables and the third is an initial guess.

时光瘦了 2025-01-01 03:36:05

fsolve 的替代方法是 < code>root

import numpy as np
from scipy.optimize import root    

def your_funcs(X):

    x, y = X
    # all RHS have to be 0
    f = [x + y**2 - 4,
         np.exp(x) + x * y - 3]

    return f

sol = root(your_funcs, [1.0, 1.0])
print(sol.x)

则会打印此信息。

[0.62034452 1.83838393]

如果您随后检查

print(your_funcs(sol.x))

[4.4508396968012676e-11, -1.0512035686360832e-11]

确认解决方案正确,

An alternative to fsolve is root:

import numpy as np
from scipy.optimize import root    

def your_funcs(X):

    x, y = X
    # all RHS have to be 0
    f = [x + y**2 - 4,
         np.exp(x) + x * y - 3]

    return f

sol = root(your_funcs, [1.0, 1.0])
print(sol.x)

This will print

[0.62034452 1.83838393]

If you then check

print(your_funcs(sol.x))

you obtain

[4.4508396968012676e-11, -1.0512035686360832e-11]

confirming that the solution is correct.

大姐,你呐 2025-01-01 03:36:05

试试这个,我向你保证它会完美地工作。

    import scipy.optimize as opt
    from numpy import exp
    import timeit

    st1 = timeit.default_timer()

    def f(variables) :
        (x,y) = variables

        first_eq = x + y**2 -4
        second_eq = exp(x) + x*y - 3
        return [first_eq, second_eq]

    solution = opt.fsolve(f, (0.1,1) )
    print(solution)


    st2 = timeit.default_timer()
    print("RUN TIME : {0}".format(st2-st1))

->

[ 0.62034452  1.83838393]
RUN TIME : 0.0009331008900937708

供参考。如上所述,您还可以通过将“fsolve”替换为“broyden1”来使用“布罗伊登近似”。有用。我做到了。

我不知道 Broyden 近似是如何工作的,但它花了 0.02 秒。

而且我建议你不要使用Sympy的功能<-确实方便,但是就速度而言,它相当慢。你会看到的。

Try this one, I assure you that it will work perfectly.

    import scipy.optimize as opt
    from numpy import exp
    import timeit

    st1 = timeit.default_timer()

    def f(variables) :
        (x,y) = variables

        first_eq = x + y**2 -4
        second_eq = exp(x) + x*y - 3
        return [first_eq, second_eq]

    solution = opt.fsolve(f, (0.1,1) )
    print(solution)


    st2 = timeit.default_timer()
    print("RUN TIME : {0}".format(st2-st1))

->

[ 0.62034452  1.83838393]
RUN TIME : 0.0009331008900937708

FYI. as mentioned above, you can also use 'Broyden's approximation' by replacing 'fsolve' with 'broyden1'. It works. I did it.

I don't know exactly how Broyden's approximation works, but it took 0.02 s.

And I recommend you do not use Sympy's functions <- convenient indeed, but in terms of speed, it's quite slow. You will see.

渔村楼浪 2025-01-01 03:36:05

我在 IDL 中使用了 Broyden 的方法来处理耦合非线性方程(通常涉及多项式和指数),但我还没有在 Python 中尝试过:

http://docs.scipy.org/doc/scipy/reference/ generated/scipy.optimize.broyden1.html#scipy.optimize.broyden1

scipy.optimize.broyden1

scipy.optimize.broyden1(F, xin, iter=None, alpha=None,duction_method='restart', max_rank=None, verbose=False, maxiter=None, f_tol=None, f_rtol=None, x_tol=无, x_rtol=无, tol_norm=无, line_search='armijo',回调=无,**kw)[来源]

使用 Broyden 的第一个雅可比近似求函数的根。

该方法也被称为“布罗伊登好方法”。

I got Broyden's method to work for coupled non-linear equations (generally involving polynomials and exponentials) in IDL, but I haven't tried it in Python:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.broyden1.html#scipy.optimize.broyden1

scipy.optimize.broyden1

scipy.optimize.broyden1(F, xin, iter=None, alpha=None, reduction_method='restart', max_rank=None, verbose=False, maxiter=None, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None, tol_norm=None, line_search='armijo', callback=None, **kw)[source]

Find a root of a function, using Broyden’s first Jacobian approximation.

This method is also known as “Broyden’s good method”.

拥抱影子 2025-01-01 03:36:05

您可以使用 openopt 包及其 NLP 方法。它有许多动态规划算法来求解非线性代数方程,包括:
goldenSection、scipy_fminbound、scipy_bfgs、scipy_cg、scipy_ncg、amsg2p、scipy_lbfgsb、scipy_tnc、bobyqa、ralg、ipopt、scipy_slsqp、scipy_cobyla、lincher、algencan可供选择。
后面的一些算法可以解决约束非线性规划问题。
因此,您可以使用如下函数将方程组引入 openopt.NLP()

lambda x: x[0] + x[1]**2 - 4、np.exp(x[0]) + x[0]*x[1]

You can use openopt package and its NLP method. It has many dynamic programming algorithms to solve nonlinear algebraic equations consisting:
goldenSection, scipy_fminbound, scipy_bfgs, scipy_cg, scipy_ncg, amsg2p, scipy_lbfgsb, scipy_tnc, bobyqa, ralg, ipopt, scipy_slsqp, scipy_cobyla, lincher, algencan, which you can choose from.
Some of the latter algorithms can solve constrained nonlinear programming problem.
So, you can introduce your system of equations to openopt.NLP() with a function like this:

lambda x: x[0] + x[1]**2 - 4, np.exp(x[0]) + x[0]*x[1]

情绪 2025-01-01 03:36:05
from scipy.optimize import fsolve

def double_solve(f1,f2,x0,y0):
    func = lambda x: [f1(x[0], x[1]), f2(x[0], x[1])]
    return fsolve(func,[x0,y0])

def n_solve(functions,variables):
    func = lambda x: [ f(*x) for f in functions]
    return fsolve(func, variables)

f1 = lambda x,y : x**2+y**2-1
f2 = lambda x,y : x-y

res = double_solve(f1,f2,1,0)
res = n_solve([f1,f2],[1.0,0.0])
from scipy.optimize import fsolve

def double_solve(f1,f2,x0,y0):
    func = lambda x: [f1(x[0], x[1]), f2(x[0], x[1])]
    return fsolve(func,[x0,y0])

def n_solve(functions,variables):
    func = lambda x: [ f(*x) for f in functions]
    return fsolve(func, variables)

f1 = lambda x,y : x**2+y**2-1
f2 = lambda x,y : x-y

res = double_solve(f1,f2,1,0)
res = n_solve([f1,f2],[1.0,0.0])
猥︴琐丶欲为 2025-01-01 03:36:05

您可以使用sympynsolve,意思是数值求解器

下面是一个使用 sympy 求解方程的示例片段:

from sympy import *

L = 4.11 * 10 ** 5
nu = 1
rho = 0.8175
mu = 2.88 * 10 ** -6
dP = 20000
eps = 4.6 * 10 ** -5

Re, D, f = symbols('Re, D, f')

nsolve((Eq(Re, rho * nu * D / mu),
       Eq(dP, f * L / D * rho * nu ** 2 / 2),
       Eq(1 / sqrt(f), -1.8 * log ( (eps / D / 3.) ** 1.11 + 6.9 / Re))),
      (Re, D, f), (1123, -1231, -1000))

其中 (1123, -1231, -1000) 是求根的初始向量。它给出:

[    13602.9381926005 - 5.34318284637546 * 10^-27 i 
   0.0479222776693451 - 1.88236901455499 * 10^-32 i
  0.00570516040320484 + 2.76709731583615 * 10^-34 i ]

虚部非常小,均为 10^(-20),因此我们可以将它们视为零,这意味着根都是实数。 Re~13602.938,D~0.047922,f~0.0057。

您可以编辑代码来求解您想要的方程,即
Eq(..., ...) 部分。不要忘记声明元素
首先使用 symbols('...')

You can use nsolve of sympy, meaning numerical solver.

Here is an example snippet which solves a equation with sympy:

from sympy import *

L = 4.11 * 10 ** 5
nu = 1
rho = 0.8175
mu = 2.88 * 10 ** -6
dP = 20000
eps = 4.6 * 10 ** -5

Re, D, f = symbols('Re, D, f')

nsolve((Eq(Re, rho * nu * D / mu),
       Eq(dP, f * L / D * rho * nu ** 2 / 2),
       Eq(1 / sqrt(f), -1.8 * log ( (eps / D / 3.) ** 1.11 + 6.9 / Re))),
      (Re, D, f), (1123, -1231, -1000))

where (1123, -1231, -1000) is the initial vector to find the root. And it gives out:

[    13602.9381926005 - 5.34318284637546 * 10^-27 i 
   0.0479222776693451 - 1.88236901455499 * 10^-32 i
  0.00570516040320484 + 2.76709731583615 * 10^-34 i ]

The imaginary part are very small, both at 10^(-20), so we can consider them zero, which means the roots are all real. Re ~ 13602.938, D ~ 0.047922 and f~0.0057.

You can edit the code to solve the equations you want, that
is the Eq(..., ...) part. Don't forget to declare elements
with symbols('...') first.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文