python - 使用牛顿法求解 Python 中的非线性方程组

Question

我正在尝试为 College 解决这个练习。我已经提交了下面的代码。但是，我对此并不完全满意。

任务是构建牛顿法的实现来求解以下非线性方程组:

为了学习牛顿法，除了类(class)之外，我还观看了这个 YouTube 视频:https://www.youtube.com/watch?v=zPDp_ewoyhM

视频中的人解释了牛顿方法背后的数学过程，并手动进行了两次迭代。

我为此做了一个 Python 实现，代码对于视频中的示例来说运行良好。尽管如此，视频中的示例涉及 2 个变量，而我的作业涉及 3 个变量。因此，我对其进行了改编。

这是代码:

import numpy as np

#### example from youtube https://www.youtube.com/watch?v=zPDp_ewoyhM

def jacobian_example(x,y):

    return [[1,2],[2*x,8*y]]

def function_example(x,y):

    return [(-1)*(x+(2*y)-2),(-1)*((x**2)+(4*(y**2))-4)]
####################################################################


### agora com os dados do exercício

def jacobian_exercise(x,y,z):

    return [[1,1,1],[2*x,2*y,2*z],[np.exp(x),x,-x]]

#print (jacobian_exercise(1,2,3))
jotinha  = (jacobian_exercise(1,2,3))

def function_exercise(x,y,z):

    return [x+y+z-3, (x**2)+(y**2)+(z**2)-5,(np.exp(x))+(x*y)-(x*z)-1]

#print (function_exercise(1,2,3))
bezao = (function_exercise(1,2,3))

def x_delta_by_gauss(J,b):

    return np.linalg.solve(J,b)

print (x_delta_by_gauss(jotinha, bezao))
x_delta_test = x_delta_by_gauss(jotinha,bezao)

def x_plus_1(x_delta,x_previous):

    x_next = x_previous + x_delta

    return x_next

print (x_plus_1(x_delta_test,[1,2,3]))

def newton_method(x_init):

    first = x_init[0]

    second = x_init[1]

    third = x_init[2]

    jacobian = jacobian_exercise(first, second, third)

    vector_b_f_output = function_exercise(first, second, third)

    x_delta = x_delta_by_gauss(jacobian, vector_b_f_output)

    x_plus_1 = x_delta + x_init

    return x_plus_1

def iterative_newton(x_init):

    counter = 0

    x_old = x_init
    print ("x_old", x_old)

    x_new = newton_method(x_old)
    print ("x_new", x_new)

    diff = np.linalg.norm(x_old-x_new)
    print (diff)

    while diff>0.000000000000000000000000000000000001:

        counter += 1

        print ("x_old", x_old)
        x_new = newton_method(x_old)
        print ("x_new", x_new)

        diff = np.linalg.norm(x_old-x_new)
        print (diff)

        x_old = x_new

    convergent_val = x_new
    print (counter)

    return convergent_val

#print (iterative_newton([1,2]))
print (iterative_newton([0,1,2]))

我很确定这段代码绝对不是完全错误的。如果我将初始值作为向量 [0,1,2] 输入，我的代码将作为输出 [0,1,2] 返回。这是一个正确的答案，它解决了上面的三个方程。

此外，如果输入 [0,2,1]，稍微不同的输入，代码也有效，它返回的答案也是正确的。

但是，如果我将初始值更改为 [1,2,3] 之类的值，我会得到一个奇怪的结果:527.7482、-1.63 和 2.14。

这个结果没有任何意义。看第一个等式，如果输入这些值，很容易看出 (527)+(-1.63)+(2.14) 不等于 3。这是错误的。

如果我将输入值更改为接近正确的解决方案，例如 [0.1,1.1,2.1]，它也会崩溃。

OK，牛顿法不能保证正确收敛。我知道。这取决于初始值等。

我的实现有任何问题吗？或者向量 [1,2,3] 只是一个“坏”初始值？

谢谢。

Answer 1

为了使您的代码更具可读性，我建议减少函数定义的数量。它们掩盖了正在发生的相对简单的计算。

我重写了我自己的版本:

def iter_newton(X,function,jacobian,imax = 1e6,tol = 1e-5):
    for i in range(int(imax)):
        J = jacobian(X) # calculate jacobian J = df(X)/dY(X) 
        Y = function(X) # calculate function Y = f(X)
        dX = np.linalg.solve(J,Y) # solve for increment from JdX = Y 
        X -= dX # step X by dX 
        if np.linalg.norm(dX)<tol: # break if converged
            print('converged.')
            break
    return X

我没有发现相同的行为:

>>>X_0 = np.array([1,2,3],dtype=float)
>>>iter_newton(X_0,function_exercise,jacobian_exercise)
converged.
array([9.26836542e-18, 2.00000000e+00, 1.00000000e+00])

甚至适用于更糟糕的猜测

>>>X_0 = np.array([13.4,-2,31],dtype=float)
>>>iter_newton(X_0,function_exercise,jacobian_exercise)
converged.
array([1.59654153e-18, 2.00000000e+00, 1.00000000e+00])