用同调矩阵求根

Question

编辑：我找到了一个可行的解决方案，但我还是希望能对这里发生的事情有更多的解释：

from scipy import optimize
from sympy import lambdify, DeferredVector

v = DeferredVector('v')
f_expr = (v[0] ** 2 + v[1] ** 2)
f = lambdify(v, f_expr, 'numpy')

zero = optimize.root(f, x0=[0, 0], method='krylov')
zero

---

原始问题：

下面是由表达式f1(x1, x2)和{}组成的矩阵M。我想知道x1和{}在M = [f1, f2] = [0, 0]时的值。在

下面的代码正在工作，除去被注释掉的根查找行。在

^{pr2}$

Answer 1

这是因为以下几个原因：

1. scipy.optimize通常坚持输入函数只包含一个参数。包装器避免了这一点，如这里所示。

2. SciPy pack中的一些线性代数函数表示输入参数与输入函数的返回值匹配。

以下代码正在工作：

import numpy as np
import pandas as pd
import sympy as sp
import matplotlib.pyplot as plt
from scipy import optimize
from sympy import init_printing, symbols, lambdify, Matrix, latex
from sympy import pi, exp, log, sqrt, sin, cos, tan, sinh, cosh, tanh
from sympy.abc import a, b, c, x, y, z, r, w
x1, x2 = symbols('x1 x2')

# Expressions
f1_expr = sin(4 * pi * x1 * x2) - 2 * x2 - x1
f2_expr = ((4 * pi  - 1) / (4 * pi)) * (exp(2 * x1) - exp(1)) + 4 * exp(1) * (x2 ** 2) - 2 * exp(1) * x1
f_expr = [f1_expr, f2_expr]

# Expressions -> NumPy function
f1 = lambdify((x1, x2), f1_expr, 'numpy')
f2 = lambdify((x1, x2), f2_expr, 'numpy')
f = np.array([f1, f2])

def _f(args):
    return [f1(args[0], args[1]), f2(args[0], args[1])]

# Matrix and it's Jacobian
M_expr = Matrix([f1_expr, f2_expr])
M_jacob_expr = M_expr.jacobian([x1, x2])

# Matrix -> NumPy function
M = lambdify((x1, x2), M_expr, [{'ImmutableMatrix': np.array}, "numpy"])
M_jacob = lambdify((x1, x2), M_jacob_expr, [{'ImmutableMatrix': np.array}, "numpy"])

# Data points
x1pts = np.arange(-2, 3, 0.01)
x2pts = np.arange(-2, 3, 0.01)
xx1pts, xx2pts = np.meshgrid(x1pts, x2pts)
# Solved over ranges for plots
z1, z2 = M(xx1pts, xx2pts)
z1 = z1.reshape(z1.shape[1], z1.shape[2])
z2 = z2.reshape(z2.shape[1], z2.shape[2])

# Find roots
results = optimize.root(_f, x0=[-0.3, 0.05], method='Krylov')
zeros = results.get('x')
# First figure
fig, ax = plt.subplots(1)

im = ax.contourf(x1pts, x2pts, z1)
ax.scatter(zeros[0], zeros[1], linewidth=5, color='k')

ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_title('f1(x1, x2)')
plt.colorbar(im)
plt.tight_layout()
plt.set_cmap('seismic')
fig.savefig('img30_1.png')
plt.show()
plt.close(fig)

# Second figure
fig, ax = plt.subplots(1)

im = ax.contourf(x1pts, x2pts, z2)
ax.scatter(zeros[0], zeros[1], linewidth=5, color='white')

ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_title('f2(x1, x2)')
plt.colorbar(im)
plt.tight_layout()
plt.set_cmap('seismic')
fig.savefig('img30_2.png')
plt.show()
plt.close(fig)