首页 文章

用python求解一个超越方程组

提问于
浏览
3

假设我有以下四个方程式:

  • cos(x)/ x = a

  • cos(y)/ y = b

  • a b = 1

  • c sinc(x)= d sinc(y)

对于未知变量 x, y, ab . 请注意 cos(x)/x=a 有多个解决方案 . 类似于变量 y . 我只对 xy 值感兴趣,这是第一个正面的根(如果重要的话) .

您可以安全地假设 a, b, cd 是已知的实常数,都是正数 .

在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}}]

因此返回

{x -> 1.31636, y -> 1.29664, a -> 0.456034, b -> 0.543966}

虽然这很容易,但我不知道如何在python中做这样的事情 . 所以,如果有人可以指导我(或只是告诉我如何)来解决这个问题,我将非常感激 .

1 回答

  • 4

    你可以使用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)
    

    它给你“旧”输出:

    [ 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]
    

相关问题