求解python中的非线性方程

Question

我有三个未知数4非线性方程X,Y和Z我想要解决.方程的形式如下:

F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ

...其中a,b并且c是常数,它们取决于F四个方程中的每个值.

解决这个问题的最佳方法是什么？

Answer 1

有两种方法可以做到这一点.

1. 使用非线性求解器
2. 线性化问题并以最小二乘意义解决问题

建立

所以,正如我理解你的问题,你知道4个不同点的F,a,b和c,你想要反演模型参数X,Y和Z.我们有3个未知数和4个观测数据点,所以问题是超定的.因此,我们将在最小二乘意义上解决.

在这种情况下使用相反的术语更常见,所以让我们翻转你的方程式.代替:

F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ

让我们写:

F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)

我们知道F,X,Y,并Z在4个不同点(例如F_0, F_1, ... F_i).

我们只是改变变量的名称,而不是方程本身.(这比我更容易思考.)

线性解决方案

实际上可以将这个等式线性化.您可以轻松地解决a^2,b^2,a b cos(c),和a b sin(c).为了使这更容易,让我们再次重新考虑:

d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)

现在这个等式简单得多:F_i = d + e X_i + f Y_i + g Z_i.这很容易做到最小二乘线性反演d,e,f,和g.然后,我们可以得到a,b以及c来自:

a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)

好的,让我们以矩阵形式写出来.我们将翻译4个观察结果(我们将编写的代码将采用任意数量的观察结果,但现在让它保持具体):

F_i = d + e X_i + f Y_i + g Z_i

成:

|F_0|   |1, X_0, Y_0, Z_0|   |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2|   |1, X_2, Y_2, Z_2|   |f|
|F_3|   |1, X_3, Y_3, Z_3|   |g|

或者:( F = G * m我是地球物理学家,所以我们使用G"格林函数"和m"模型参数".通常我们也用d"数据"代替F.)

在python中,这将转换为:

def invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

非线性解决方案

scipy.optimize正如@Joe建议的那样,您也可以使用它来解决这个问题.最容易使用的函数scipy.optimize是scipy.optimize.curve_fit默认使用Levenberg-Marquardt方法.

Levenberg-Marquardt是一种"爬山"算法(在这种情况下,它会下坡,但无论如何都会使用该术语).从某种意义上说,您可以对模型参数进行初步猜测(默认情况下为全部scipy.optimize),然后按照observed - predicted参数空间的斜率下坡到底部.

警告:采用正确的非线性反演方法,初始猜测和调整方法的参数是非常"黑暗的艺术".你只能通过这样做来学习它,并且在很多情况下,事情将无法正常工作.如果你的参数空间相当平滑(这应该是),Levenberg-Marquardt是一个很好的通用方法.除了更常见的方法(如模拟退火)之外,还有很多其他方法(包括遗传算法,神经网络等)在其他情况下更好.我不打算在这里深入研究这一部分.

有一个常见的问题,一些优化工具包试图纠正,而scipy.optimize不是试图处理.如果您的模型参数具有不同的大小(例如a=1, b=1000, c=1e-8),则需要重新调整尺寸以使它们的大小相似.否则scipy.optimize,"爬山"算法(如LM)将无法准确计算出局部梯度的估计值,并且会产生非常不准确的结果.现在,我假设a,b以及c具有相对类似的量值.此外,请注意,基本上所有非线性方法都需要您进行初始猜测,并且对该猜测很敏感.我将它留在下面(只是将其作为p0kwarg传递给curve_fit),因为默认值a, b, c = 1, 1, 1是相当准确的猜测a, b, c = 3, 2, 1.

随着注意事项的出现,curve_fit期望传递一个函数,一组观察点(作为单个ndim x npoints数组)和观察值.

所以,如果我们写这样的函数:

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

我们需要将其包装以接受稍微不同的参数,然后再传递给它curve_fit.

简而言之:

def nonlinear_invert(f, x, y, z):
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

这两种方法的独立示例:

为了给你一个完整的实现,这是一个例子

1. 生成随机分布的点来评估函数,
2. 评估这些点上的函数(使用set model parameters),
3. 为结果增加噪音,
4. 然后使用上述线性和非线性方法反演模型参数.

import numpy as np
import scipy.optimize as opt

def main():
    nobservations = 4
    a, b, c = 3.0, 2.0, 1.0
    f, x, y, z = generate_data(nobservations, a, b, c)

    print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
    print linear_invert(f, x, y, z)

    print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
    print nonlinear_invert(f, x, y, z)

def generate_data(nobservations, a, b, c, noise_level=0.01):
    x, y, z = np.random.random((3, nobservations))
    noise = noise_level * np.random.normal(0, noise_level, nobservations)
    f = func(x, y, z, a, b, c) + noise
    return f, x, y, z

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

def linear_invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

def nonlinear_invert(f, x, y, z):
    # "curve_fit" expects the function to take a slightly different form...
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

main()