python - 如何绘制梯度(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,是)；这是代码的样子，但它报告错误

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 符号与 float 相乘。

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

更新

@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) = Sum over n and m of [Cn * Dm * sin(2pinx) * cos(2pimy)]; sum goes from -10 to 10; Cn and Dm are coefficients, and i calculated that CkDl = sin(2pik)/(k^2 +l^2) (i used here k and l as one of the indices from the sum over n and 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), with
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 的定义中也可以看出)您的函数是周期性的，因此如果您以更小的限制和更少的网格绘制相同的东西点，您会看到更多特征和更大的箭头。