SymPy: создание функции numpy из диагональной матрицы, которая принимает массив numpy

Основываясь на примере, который я нашел здесь, я пытаюсь создать функцию из диагональной матрицы, которая была создана с помощью sumpy.diag

myM = Matrix([
[x1, 4, 4],
[4, x2, 4],
[4, 4, x3]])

Если это было создано с использованием этой процедуры, например:

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
x3 = sp.Symbol('x3')
X = sp.Matrix([x1, x2, x3])

myM = 4 * sp.ones(3, 3)
sp.diag(*X) + myM - sp.diag(*np.diag(myM))

теперь я хотел бы создать функцию, используя lambdify of ufuncify, которая принимает в качестве ввода значение numpy.array или length 3 (например, np.array([0.1,0.2,0.3])) и дает результат в виде матрицы в соответствии с myM

myM = Matrix([
[0.1, 4, 4],
[4, 0.2, 4],
[4, 4, 0.3]])

В конце концов мне нужно создать матрицу якобиана, символически используя этот метод: И так как функциональная форма может меняться во время вычисления, то использование Якобиана, рассчитанное символически, было бы очень полезно.

Ответ 1

Создание числовой матрицы 3 на 3 из числового вектора на самом деле не является символом SymPy, поскольку никакие символы не задействованы. Рассмотрим следующее, где аргумент d представляет собой массив, содержащий диагональные элементы.

def mat(d):
    return np.diag(d-4) + 4

Вышеуказанная функция возвращает массив 2d NumPy. Чтобы вернуть матрицу SymPy, используйте

def mat(d):
    return sp.Matrix(np.diag(d-4) + 4)

Когда d имеет чрезвычайно малые значения, вычитание с последующим добавлением может привести к потере точности: например, (1e-20 - 4) + 4 оценивается до нуля. Более безопасная альтернатива -

def mat(d):
    diagmat = np.diag(d) 
    return diagmat + np.fromfunction(lambda i, j: (i != j)*4, diagmat.shape)

Ответ 2

вы можете .subs() плавать значения в соответствующие символы:

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
x3 = sp.Symbol('x3')
X = sp.Matrix([x1, x2, x3])

myM = 4 * sp.ones(3, 3)
smyM=sp.diag(*X) + myM - sp.diag(*np.diag(myM))

fcoefs = [(a, f) for a, f in (zip([x1, x2, x3], np.array([0.1,0.2,0.3])))]

fmyM = smyM.subs(fcoefs)

smyM
Out[105]: 
Matrix([
[x1,  4,  4],
[ 4, x2,  4],
[ 4,  4, x3]])

fmyM
Out[106]: 
Matrix([
[0.1,   4,   4],
[  4, 0.2,   4],
[  4,   4, 0.3]])

кажется прекрасной матрицей sympy.matrices.dense.MutableDenseMatrix после:

fmyM @ myM
Out[107]: 
Matrix([
[32.4, 32.4, 32.4],
[32.8, 32.8, 32.8],
[33.2, 33.2, 33.2]])

может потребоваться преобразование в np.array для полного использования с numpy

ниже приведена часть моего кода, показывающая больше шаблона, который я использовал:

def ysolv(coeffs):
    x,y,a,b,c,d,e = symbols('x y a b c d e')
    ellipse = a*y**2 + b*x*y + c*x + d*y + e - x**2
    y_sols = solve(ellipse, y)
    print(*y_sols, sep='\n')

    num_coefs = [(a, f) for a, f in (zip([a,b,c,d,e], coeffs))]
    y_solsf0 = y_sols[0].subs(num_coefs)
    y_solsf1 = y_sols[1].subs(num_coefs)

    f0 = lambdify([x], y_solsf0)
    f1 = lambdify([x], y_solsf1)
    return f0, f1

f0, f1 = ysolv(t[0])

y0 = [f0(x) for x in xs]
y1 = [f1(x) for x in xs]
...

from: fooobar.com/info/504421/... (да, мой "feeloutXrange" есть так плохо, что он должен быть показан)