2

我最近刚刚学会了如何使用 SciPy 的四元积分来积分一个二变量函数，仅针对一个变量。例如，假设我们有一个函数 f(x,y)。我对 x 进行积分，以获得 f(y)。然后我绘制了 f(y) 与 y 值的 NumPy 数组。这里没问题。但是，我现在想将我的原始被积函数乘以雅可比矩阵（雅可比矩阵的各个条目中的每个函数也是 x 和 y 的函数）。

我现在无法完全概念化的问题是如何将部分微分合并到我的原始代码中——它非常干净地集成了我最初的二变量函数。我知道 SymPy 是符号积分的一个很好的选择，但我无法弄清楚如何评估雅可比行列式（因为这需要导数）以一种与其结果相乘的方式相兼容的方式乘以一个单独的函数，然后集成该产品。是否有一种形式的lambdify，它不是将我的 SymPy 函数转换为它们的 NumPy 对应函数，而是将 SymPy 函数转换为可用于 SciPy 的四元集成的函数？还是他们的另一种常见做法？

from sympy import symbols, diff
import numpy as np
from scipy.integrate import quad
from scipy import integrate
import matplotlib.pyplot as plt

## define f(x,y) function *** BUT INTEGRATING dx ***
## where the bounds are as a function of y ***
## SciPy and NumPy are used in these steps

def example(y):
    return lambda x: (x * y)

def example_integral(y):
    return quad(example(y), 0, y, args=())


## F(y) = (y^3) / 2

print('F(y) --> F(4) = (4^3)/2 = ',example_integral(4)[0])
print('Integration works')


## Plotting f(y) vs. y
## This graph should follow cubic power
example_range = np.linspace(-11,10,20)
example_array = []

for i in example_range:
    example_array.append(example_integral(i)[0])

plt.figure(figsize=(6,6))    
plt.plot(example_range,example_array)
plt.title("Example Plot of original integrated function\nF(y) vs y")
plt.show()

## Here I will differentiate two separate functions for the Jacobian
## SymPy is used in these following steps
x, y = symbols('x y', real=True)
f = (2*x*(y**2)) + y**3
g = x/y

df_dx = diff(f, x)
df_dy = diff(f, y)
dg_dx = diff(g, x)
dg_dy = diff(g, y)

## I compute the Jacobian here
jacobian = (df_dy * dg_dx) - (df_dx * dg_dy)
print("The resulting Jacobian is: ", jacobian)

根据建议，我包含了最能解释我的困境的示例代码。如果这不遵守“最小示例”指南，请发表评论，我会解决这个问题。如您所见，在第一个块中，我使用 SciPy 和 NumPy 来集成和绘制我的原始函数。在第二个块中，我使用 SymPy 来区分两个独立的函数并计算雅可比行列式。最终，我想将得到的雅可比行列式（在第二个块中）乘以我的原始函数（在第一个块中）。在我有了这个新功能（两者的产物）之后，我想简单地执行相同的集成过程。附在这里 --> 1是我尝试运行此想法时收到的错误消息的屏幕截图（这是预期的）。如果我能做些什么来更好地解释这一点，请告诉我。谢谢！

4

1 回答 1

1

我只能大胆猜测 SymPyufuncify会帮助你。它将 SymPy 表达式转换为向量化的 NumPy ufuncs，可以以向量化的方式应用。因此它也可以被 SciPy 的正交例程使用。

于 2020-05-11T06:49:00.670 回答