可以使用scipy.optimize的Newton-Raphson方法来处理多变量系统吗?

huangapple go评论64阅读模式
英文:

Can i use the Newton-Raphson method of scipy.optimize with a multiple variable system?

问题

我想用 scipy 进行多变量系统的 Newton-Raphson 方法。所以我按照这个文档。这是我尝试解决问题的示例代码:

import numpy as np
from scipy.optimize import newton

def f(x):
    F = [0, 0]
    F[0] = 5 * x[0] + 7 * x[1] - 6
    F[1] = 10 * x[0] - 3 * x[1] - 46
    return F

def jacobian(x):
    F = np.zeros((2, 2))
    F[0, 0] = 5
    F[0, 1] = 7
    F[1, 0] = 10
    F[1, 1] = -3
    return F

q = [10, -10]
result = newton(f, q, fprime=jacobian, maxiter=10000, full_output=1)
print(result)

当我在不使用 'fprime' 参数的情况下运行此程序时,我得到了正确的答案:

result(root=array([ 3.99999999, -2.00000004]), converged=array([ True, True]), zero_der=array([False, False]))

但这是割线法,对于更复杂的问题,我需要使用牛顿法。所以当我使用 'fprime' 参数运行代码时,我得到了以下错误:

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

看到这个错误,我尝试了几种解决方案,比如将 jacobian 函数中的矩阵 F 更改为普通列表:

def jacobian(x):
    F = [0, 0, 0, 0]
    F[0] = 5
    F[1] = 7
    F[2] = 10
    F[3] = -3
    return F

或将所有变量更改为 numpy 数组,但我仍然得到类似的索引错误。

我想知道是否可以使用这个 scipy 函数解决这个方程组,如果不行,我想了解其他在 Python 中的替代方案。

英文:

I wanted to do the Newton-Raphson method with scipy with a multivariable system. So, i followed this documentation. And here is an example code where i tried to solved my problem:

import numpy as np
from scipy.optimize import newton


def f(x):
    F = [0, 0]
    F[0] = 5 * x[0] + 7 * x[1] - 6
    F[1] = 10 * x[0] - 3 * x[1] - 46

    return F


def jacobian(x):
    F = np.zeros((2, 2))
    F[0, 0] = 5
    F[0, 1] = 7
    F[1, 0] = 10
    F[1, 1] = -3

    return F


q = [10, -10]
result = newton(f, q, fprime=jacobian, maxiter=10000, full_output=1)
print(result)

When i run this program without the 'fprime' parameter (in the newton function), i get the correct answer:

> result(root=array([ 3.99999999, -2.00000004]), converged=array([ True, True]), zero_der=array([False, False]))

But this is the secant method and i need the Newton method for a more complex problem. So when i run the code with the parameter 'fprime' i get the following error:

> IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

Seeing this error i tried several solutions, like changing the matrix F from the jacobian function to a normal list:

def jacobian(x):
    F = [0,0,0,0]
    F[0] = 5
    F[1] = 7
    F[2] = 10
    F[3] = -3

    return F

Or changing all the variable to numpy arrays, but i keep getting similar errors with index.

I want know if it is posible to do this system of equations with this scipy function and if not, i would like other alternatives with python.

答案1

得分: 0

The traceback is

  Cell In[46], line 26
    result = newton(f, q, fprime=jacobian, maxiter=10000, full_output=1)

  File ~\miniconda3\lib\site-packages\scipy\optimize\_zeros_py.py:277 in newton
    return _array_newton(func, x0, fprime, args, tol, maxiter, fprime2,

  File ~\miniconda3\lib\site-packages\scipy\optimize\_zeros_py.py:401 in _array_newton
    dp = fval[nz_der] / fder[nz_der]

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

Changing the ipython %xmode` to verbose, the error line, with context and variables is:

File ~\miniconda3\lib\site-packages\scipy\optimize\_zeros_py.py:401, in _array_newton(func=<function f>, x0=[10, -10], fprime=<function jacobian>, args=(), tol=1.48e-08, maxiter=10000, fprime2=None, full_output=1)
    399     break
    400 # Newton step
    400 # Newton step
--> 401 dp = fval[nz_der] / fder[nz_der]
        nz_der = array([[ True,  True],
       [ True,  True]])
        fval = array([-26,  84])
        fder = array([[ 5.,  7.],
       [10., -3.]])

I think fval is the 1d array produced by f(q). fder is the 2d array produced by jacobian, and nz_der looks like a matching shape boolean.

You/we need to review the newton docs to see what shape of a result the fprime must be.

The docs aren't specific about what fprime must return, but from that error, and the examples like fprime=lambda x: 3 * x**2), it must be a 1d array like f(q), the same shape as in guess q. You are trying to provide the 2d dz/dx.

Using np.array([5,7,10,-3]) has indexing problem at the same place, only with a (4,) boolean instead of the (2,2).

From your linked docs:

> If x0 is a sequence with more than one item, newton returns an array: the zeros of the function from each (scalar) starting point in x0. In this case, func must be vectorized to return a sequence or array of the same shape as its first argument. If fprime (fprime2) is given, then its return must also have the same shape: each element is the first (second) derivative of func with respect to its only variable evaluated at each element of its first argument.

英文:

The traceback is

  Cell In[46], line 26
    result = newton(f, q, fprime=jacobian, maxiter=10000, full_output=1)

  File ~\miniconda3\lib\site-packages\scipy\optimize\_zeros_py.py:277 in newton
    return _array_newton(func, x0, fprime, args, tol, maxiter, fprime2,

  File ~\miniconda3\lib\site-packages\scipy\optimize\_zeros_py.py:401 in _array_newton
    dp = fval[nz_der] / fder[nz_der]

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

Changing the ipython %xmode` to verbose, the error line, with context and variables is:

File ~\miniconda3\lib\site-packages\scipy\optimize\_zeros_py.py:401, in _array_newton(func=<function f>, x0=[10, -10], fprime=<function jacobian>, args=(), tol=1.48e-08, maxiter=10000, fprime2=None, full_output=1)
    399     break
    400 # Newton step
    400 # Newton step
--> 401 dp = fval[nz_der] / fder[nz_der]
        nz_der = array([[ True,  True],
       [ True,  True]])
        fval = array([-26,  84])
        fder = array([[ 5.,  7.],
       [10., -3.]])

I think fval is the 1d array produced by f(q). fder is the 2d array produced by jacobian, and nz_der looks like a matching shape boolean.

You/we need to review the newton docs to see what shape of a result the fprime must be.

The docs aren't specific about what fprime must return, but from that error, and the examples like fprime=lambda x: 3 * x**2), it must be a 1d array like f(q), the same shape as in guess q. You are trying to provide the 2d dz/dx.

Using np.array([5,7,10,-3]) has indexing problem at the same place, only with a (4,) boolean instead of the (2,2).

From your linked docs:

> If x0 is a sequence with more than one item, newton returns an array: the zeros of the function from each (scalar) starting point in x0. In this case, func must be vectorized to return a sequence or array of the same shape as its first argument. If fprime (fprime2) is given, then its return must also have the same shape: each element is the first (second) derivative of func with respect to its only variable evaluated at each element of its first argument.

huangapple
  • 本文由 发表于 2023年3月31日 22:38:51
  • 转载请务必保留本文链接:https://go.coder-hub.com/75899797.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定