Численно приближенный якобиан вектора

Question

Я должен найти якобиан функции

def f(x):
    return np.array([
        np.sin(x[0])  + 0.5 * (x[0] - x[1])**3 - 1.0,
        0.5 * (x[1] - x[0])**3 + x[1]
    ])

Одно из требований: если dx не указан, мне нужно установить его на 10e-5. Мне дают два параметра для распечатки, которые будут запускать функцию, которую я пишу. Вот что у меня есть:

def jac(f, x0, dx=None):
    x0 = np.atleast_1d(x0)
    x = np.copy(x0)
    f0 = np.copy(f)
    N = x.size
    i = 0
    if dx is None:
        dx = 10^-5*np.ones(N)
    else:
       dx = np.atleast_1d(N)
    for i in range(N):
      J=np.zeros([N,N])
      x[i]=x[i]+dx[i]
      J[:,i]=(f(x)-f0)/dx[i]
    return J

print(jac(f,(1,2),dx=(1e-3,1e-3)))
print('\n')
print(jac(f, np.array((1.0,2.0))))

Я продолжаю получать сообщения об ошибках неподдерживаемых типов операндов для -: 'float' и 'function' для этого раздела кодирования

TypeError                                 Traceback (most recent call last)
<ipython-input-58-3e78218d0309> in <module>()
     23       J[:,i]=(f(x)-f0)/dx[i]
     24     return J
---> 25 print(jac(f,(1,2),dx=(1e-3,1e-3)))
     26 print('\n')
     27 #print(jac(f, np.array((1.0,2.0))))

<ipython-input-58-3e78218d0309> in jac(f, x0, dx)
     21       J=np.zeros([N,N])
     22       x[i]=x[i]+dx[i]
---> 23       J[:,i]=(f(x)-f0)/dx[i]
     24     return J
     25 print(jac(f,(1,2),dx=(1e-3,1e-3)))

TypeError: unsupported operand type(s) for -: 'float' and 'function'

Я не уверен, как решить эту ошибку или что она говорит. Буду признателен за любую помощь в исправлении ошибки или мой код в целом.

Answer 1

В вашем коде есть несколько ошибок:

• f0 - это функция, а не массив np.
• dx=np.atleast_1d(N) не имеет смысламне.
• Вы устанавливаете J в каждой итерации в качестве нулевой матрицы.

Вот способ, которым вы могли бы сделать это:

import numpy as np

f = lambda x: np.array([np.sin(x[0])  + 0.5 * (x[0] - x[1])**3 - 1.0, 0.5 * (x[1] - x[0])**3 + x[1]])

def jac(f,x0,dx=None):
    x0 = np.atleast_1d(x0)
    N  = x0.size
    J = np.zeros((N, N))
    if dx is None:
        dx=10^-5*np.ones(N)
    else:
        dx=np.atleast_1d(dx)
    for i, ei in enumerate(np.eye(N)):
        J[:,i] = (f(x0 + dx[i]*ei)-f(x0))/dx[i]
    return J

Кстативместо этого вы можете использовать недокументированную функцию ок_приемлем из scipy.optimize._numdiff. Я буду использовать его для проверки нашей реализации:

from scipy.optimize._numdiff import approx_derivative

print(jac(f,(1,2),dx=(1e-6,1e-6)))
print(approx_derivative(f,(1,2)))

дает мне

[[ 2.04030038 -1.5000015 ]
 [-1.4999985   2.5000015 ]]

[[ 2.04030231 -1.5       ]
 [-1.5         2.5       ]]