如何绘制 grad(f(x,y))？

Question

我想计算和绘制两个变量的任何标量函数的梯度。如果你真的想要一个具体的例子，让我们说 f=x^2+y^2 其中 x 从 -10 到 10 和 y 相同。如何计算和绘制 grad(f)？解决方案应该是矢量，我应该看到矢量线。我是python新手，所以请用简单的词。

编辑：

@Andras Deak：谢谢你的帖子，我尝试了你的建议，而不是你的测试函数 (fun=3*x^2-5*y^2) 我使用了我定义为 V(x, y);代码是这样的，但是报错了

import numpy as np
import math 
import sympy
import matplotlib.pyplot as plt

def V(x,y):
    t=[]
    for k in range (1,3): 
        for l in range (1,3):
            t.append(0.000001*np.sin(2*math.pi*k*0.5)/((4*(math.pi)**2)* (k**2+l**2)))  
            term = t* np.sin(2 * math.pi * k * x/0.004) * np.cos(2 * math.pi * l * y/0.004)
            return term 
    return term.sum()

x,y=sympy.symbols('x y')
fun=V(x,y)
gradfun=[sympy.diff(fun,var) for var in (x,y)]
numgradfun=sympy.lambdify([x,y],gradfun)

X,Y=np.meshgrid(np.arange(-10,11),np.arange(-10,11))
graddat=numgradfun(X,Y)
plt.figure()
plt.quiver(X,Y,graddat[0],graddat[1])
plt.show()

AttributeError: 'Mul' object has no attribute 'sin'

假设我删除了罪，我得到另一个错误：

TypeError: can't multiply sequence by non-int of type 'Mul'

我阅读了 sympy 教程，它说“符号计算系统（如 SymPy）的真正强大之处在于能够以符号方式进行各种计算”。我明白了，我只是不明白为什么我不能将 x 和 y 符号与浮点数相乘。

解决这个问题的方法是什么？ :(请帮忙！

更新

@Andras Deak：我想让事情变得更短，所以我从 V(x,y) 和 Cn*Dm 的原始公式中删除了许多常量。正如您所指出的，这导致 sin 函数始终返回 0（我刚刚注意到）。对此表示歉意。今天晚些时候，当我详细阅读您的评论时，我将更新该帖子。非常感谢！

更新 2 我更改了电压表达式中的系数，结果如下：

它看起来不错，只是箭头指向相反的方向（它们应该从红点出来并进入蓝色点）。你知道我怎么能改变它吗？如果可能的话，你能告诉我增加箭头大小的方法吗？我尝试了另一个主题 (Computing and drawing vector fields) 中的建议：

skip = (slice(None, None, 3), slice(None, None, 3))

这只会绘制每三个箭头，matplotlib 会自动缩放，但它对我不起作用（当我添加这个时，我输入的任何数字都不会发生任何事情） 你已经帮了我很大的忙，我感激不尽！

Answer 1

这是使用sympy 和numpy 的解决方案。这是我第一次使用 sympy，所以其他人会/可能会想出更好、更优雅的解决方案。

import sympy

#define symbolic vars, function
x,y=sympy.symbols('x y')
fun=3*x**2-5*y**2

#take the gradient symbolically
gradfun=[sympy.diff(fun,var) for var in (x,y)]

#turn into a bivariate lambda for numpy
numgradfun=sympy.lambdify([x,y],gradfun)

现在您可以使用numgradfun(1,3) 计算(x,y)==(1,3) 处的梯度。然后可以使用此功能进行绘图，您说可以做到。

对于绘图，您可以使用例如matplotlib 的quiver，如下所示：

import numpy as np
import matplotlib.pyplot as plt

X,Y=np.meshgrid(np.arange(-10,11),np.arange(-10,11))
graddat=numgradfun(X,Y)

plt.figure()
plt.quiver(X,Y,graddat[0],graddat[1])
plt.show()

更新

您为要计算的函数添加了规范。它包含取决于x 和y 的术语乘积，这似乎打破了我的上述解决方案。我设法想出了一个新的来满足您的需求。但是，您的功能似乎没有什么意义。从您编辑的问题：

t.append(0.000001*np.sin(2*math.pi*k*0.5)/((4*(math.pi)**2)* (k**2+l**2)))
term = t* np.sin(2 * math.pi * k * x/0.004) * np.cos(2 * math.pi * l * y/0.004)

另一方面，从您对这个答案的相应评论：

V(x,y) = [Cn * Dm * sin(2pinx) * cos(2pimy)] 的 n 和 m 之和；总和从 -10 到 10； Cn 和 Dm 是系数，我计算了 即 CkDl = sin(2pik)/(k^2 +l^2) （我在这里使用 k 和 l 作为其中之一 n 和 m 之和的索引）。

我对此有几个问题：sin(2*pi*k) 和 sin(2*pi*k/2)（对于整数 k，前因子中的两个竞争版本始终为零，在每个 (x,y) 处为您提供一个常量零 V。此外，在您的代码中，您的三角函数中有神奇的频率因子，而注释中缺少这些因子。如果您将 x 乘以 4e-3，您会大大改变函数的空间依赖性（通过将波长改变大约一千倍）。所以你应该真正决定你的功能是什么。

所以这是我假设的解决方案

V(x,y)=sum_{k,l = 1 to 10} C_{k,l} * sin(2*pi*k*x)*cos(2*pi*l*y)，有
C_{k,l}=sin(2*pi*k/4)/((4*pi^2)*(k^2+l^2))*1e-6

这是您各种版本的函数的组合，在前置因子中修改了sin(2*pi*k/4)，以便具有非零函数。我希望您在找出正确的数学模型后，能够根据您的实际需要来修正数值因子。

完整的代码如下：

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

def CD(k,l):
    #return sp.sin(2*sp.pi*k/2)/((4*sp.pi**2)*(k**2+l**2))*1e-6
    return sp.sin(2*sp.pi*k/4)/((4*sp.pi**2)*(k**2+l**2))*1e-6

def Vkl(x,y,k,l):
    return CD(k,l)*sp.sin(2*sp.pi*k*x)*sp.cos(2*sp.pi*l*y)

def V(x,y,kmax,lmax):
    k,l=sp.symbols('k l',integers=True)
    return sp.summation(Vkl(x,y,k,l),(k,1,kmax),(l,1,lmax))


#define symbolic vars, function
kmax=10
lmax=10
x,y=sp.symbols('x y')
fun=V(x,y,kmax,lmax)

#take the gradient symbolically
gradfun=[sp.diff(fun,var) for var in (x,y)]

#turn into bivariate lambda for numpy
numgradfun=sp.lambdify([x,y],gradfun,'numpy')
numfun=sp.lambdify([x,y],fun,'numpy')

#plot
X,Y=np.meshgrid(np.linspace(-10,10,51),np.linspace(-10,10,51))
graddat=numgradfun(X,Y)
fundat=numfun(X,Y)

hf=plt.figure()
hc=plt.contourf(X,Y,fundat,np.linspace(fundat.min(),fundat.max(),25))
plt.quiver(X,Y,graddat[0],graddat[1])
plt.colorbar(hc)
plt.show()

我使用一些辅助函数定义了您的 V(x,y) 函数以实现透明度。我将求和截止值保留为文字参数，kmax 和 lmax：在您的代码中它们是 3，在您的评论中它们被称为 10，无论如何它们应该是无穷大。

渐变采用与以前相同的方式，但是当使用lambdify 转换为numpy 函数时，您必须设置一个额外的字符串参数'numpy'。这将使生成的 numpy lambda 接受数组输入（本质上它将使用np.sin 而不是math.sin，cos 也是如此）。

我还将网格的定义从array 更改为np.linspace：这通常更方便。由于您的函数在整数网格点处几乎是恒定的，因此我创建了一个更密集的网格进行绘图（51 个点，同时保持 (-10,10) 的原始限制固定）。

为清楚起见，我添加了更多图表：contourf 显示函数的值（等高线应始终与梯度向量正交），以及指示函数值的颜色条。结果如下：

显然构图不是最好的，但我不想偏离你的规格太多。该图中的箭头实际上几乎不可见，但正如您所看到的（从 V 的定义中也可以看出）您的函数是周期性的，因此如果您以较小的限制和较少的网格点绘制相同的东西，您将查看更多功能和更大的箭头。