SymPy не оценивает экспоненциальную функцию с помощью subs

#python #sympy

#python #sympy

Вопрос:

Я пытаюсь написать реализацию многомерного метода Ньютона-Рафсона, в Python . Для начала я пытаюсь решить систему:

 10 * x   3 * y * y - 3
x * x - exp(y) - 2
 

Хотя я намерен распространить это на любую произвольную m x n систему (при условии, что решение возможно / существует, конечно).

Мой код выглядит следующим образом:

 from dataclasses import dataclass
from sympy import *

x, y, z = symbols('x y z')


@dataclass
class Solve:

    @staticmethod
    def newton_raphson():

        F = Matrix([10 * x   3 * y * y - 3, x * x - exp(y) - 2])
        v = Matrix([x, y])
        print(J := F.jacobian(v))
        print(J * -1)

        xx, yy = 0, 0
        for i in range(10):
            A = J.subs({"x": xx, "y": yy})
            b = F.subs({"x": xx, "y": yy})

            update = linsolve((A, b), [x, y])

            (xx, yy) = tuple(*update) # should be adding, not setting equal to

            print(xx, yy)




mySolver = Solve
mySolver.newton_raphson()
 

Моя проблема в том, что, когда я выполняю .subs() работу с матрицами F и J , вычисляется все, кроме экспоненциальной функции. То есть строка print(xx, yy) выводит:

 3*(573 - 50*exp(3))/(20*(27 - 25*exp(3))) 5*(13   20*exp(3))/(4*(-27   25*exp(3)))

(-28650000*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   9) - 7500000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   9) - 53550000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   6) - 105598890*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6)))   92498679*exp(5*(351   215*exp(3))/(4*(-1350*exp(3)   729   625*exp(6))))   164313360*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   3)   15235500*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   6)   164358750*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   3))/(-81000000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   6) - 107010000*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   6) - 31492800*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6)))   72404280*exp(5*(351   215*exp(3))/(4*(-1350*exp(3)   729   625*exp(6))))   38032200*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   3)   87480000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   3)   9000000*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   9)   25000000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   9)) (-16200000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   6) - 40086000*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   6) - 68751693*exp(5*(351   215*exp(3))/(4*(-1350*exp(3)   729   625*exp(6)))) - 6298560*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6)))   17496000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   3)   68874390*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   3)   10450000*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   9)   5000000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   9))/(8*(-2025000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   6) - 2675250*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   6) - 787320*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6)))   1810107*exp(5*(351   215*exp(3))/(4*(-1350*exp(3)   729   625*exp(6))))   950805*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   3)   2187000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   3)   225000*exp(1755/(4*(-1350*exp(3)   729   625*exp(6)))   1075*exp(3)/(4*(-1350*exp(3)   729   625*exp(6)))   9)   625000*exp(625*exp(6)/(-1350*exp(3)   729   625*exp(6))   9)))

...etc

 

Хотя мне нужно строго числовое значение, для обновления догадок xx и yy .

Комментарии:

1. В Python есть синтаксис для полномочий, поэтому вместо 3*y*y вы можете написать 3*y**2 .

Ответ №1:

Sympy специализируется на символьной оценке. exp(3) является иррациональным числом и, следовательно, не может быть представлено численно. Sympy знает это.

Чтобы принудительно выполнить численное приближение, вы можете использовать N функцию из sympy, как в

A = N(J.subs({"x": xx, "y": yy}))

Комментарии:

1. Вы также можете использовать J.evalf(subs={"x", ...}) , который потенциально более точен в числовом отношении. Обычно что-то подобное было бы сделано lambdify для ускорения (например, это nsolve то, что делает sympy).