时间戳

首页 > 解决方案 > 使用python数值求解有限定义区间上的非线性方程组

python - 使用python数值求解有限定义区间上的非线性方程组

问题描述

假设以下问题：我必须要函数，其中一个由 t 参数化，并且两者都定义在自变量 z 上：f(t, z, *ags), g(z, *args)。

我想找到一对值 (t0, z0) 以便

函数及其关于 z 的导数是相同的 f(t0, z0, *args) = g(z0, *args) 和 df/dz(t0, z0, *args) = dg/dz(z0, *args ）。

我知道存在一个解决方案，并且我得到了一个合理的起点（tS，zS）。但是，至少其中一个函数仅在指定的间隔 [zL .. zH] 上定义（我知道）。

我现在的问题是，这是在 python 中以数值方式求解方程组的最佳方法。

我尝试了 scipys fsolve，但它似乎失败了，我认为是因为它无法处理有限的定义间隔。我尝试了差分进化包，只是为了最小化复合函数，但这似乎完全是矫枉过正。

我有所有函数及其衍生物的表达式（尽管它们很复杂）。

肯定有一个简单的python rootfinder，它能够解决两个非线性方程组，它们只在有限的区间上定义？！

出于某种原因，我发现要么能够求解方程组，不考虑限制，要么能够限制但一次只能求解一个方程......

如果有人能指出我正确的方向，我将不胜感激！

标签： pythonsolvernonlinear-functions

解决方案

这个问题需要假设 1) 存在多个解决方案和 2) 至少一个解决方案在指定区间的范围内。如果不满足这些假设，您的问题就变成了优化——最小化您的函数之间的差异。

我生成了以下满足假设的示例：

import numpy as np
from scipy.optimize import fsolve

def f(t,z):
    return t**2 + z + z**2 + np.sin(z)
def g(z):
    return z**2
def dfdz(t,z):
    return 1 + 2*z + np.cos(z)
def dgdz(z):
    return 2*z

def solve(x):
    t,z = x

    #residual array
    r = np.zeros(2)
    #match equations
    r[0] = (f(t,z) - g(z))**2
    #match derivatives
    r[1] = dfdz(t,z) - dgdz(z)

    return r

x0 = [10,-10]
x = fsolve(solve,x0)
t,z = x
print(t,z)

在这些初始条件下，答案是 [3.06998012761 -9.42477800292]。通过在超出范围时添加残差，使用 fsolve 可以要求有界变量。这是一个非常 hacky 的解决方案，不是很健壮。我们可以通过以下修改将 z 绑定到 [-5,0] 的范围（这里还有另一个正确答案）：

def solve(x,zL,zH):
    t,z = x

    #penalize deviation from range
    if z < zL:
        e = z - zL
    elif z > zH:
        e = z - zH
    else:
        e = 0

    #residual array
    r = np.zeros(2)
    #match equations
    r[0] = (f(t,z) - g(z))**2 + (e)**2
    #match derivatives
    r[1] = dfdz(t,z) - dgdz(z)

    return r


x0 = [10,-10]
x = fsolve(solve,x0,args=(-5,0))
t,z = x
print(t,z)

新答案是 [1.77245385 -3.14159263]。

如前所述，有界变量在优化问题中更常见，因此优化包会更直观地处理它。在 python 中，有标准的 scipy.optimize.minimize，或者更强大的包，例如Pyomo或GEKKO。这是 GEKKO 中的等效问题：

from gekko import GEKKO

#initialize a GEKKO model
m = GEKKO()

#add GEKKO variables
t = m.Var(value = -10)
z = m.Var(value = -10, lb=-5, ub=0)

#define the constraints
m.Equation(t**2 + z + z**2 + m.sin(z) == z**2)
m.Equation(1 + 2*z + m.cos(z) == 2*z)

#solve a system of non-dynamic equations
m.options.IMODE = 1 
m.solve(disp=False)

print(t.value,z.value)

[-1.772454] [-3.142608]

推荐阅读