如何绘制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符号乘以浮点数.

这是怎么回事？:( 请帮忙!

UPDATE

@Andras Deak:我想让事情更短,所以我从V(x,y)和Cn*Dm的原始公式中删除了许多常量.正如你所指出的,这导致sin函数总是返回0(我刚注意到).为此道歉.我将在今天晚些时候更新帖子,当时我会详细阅读您的评论.十分感谢!

更新2 我在表达式中更改了电压系数,这是结果:

它看起来不错,只是箭头指向相反的方向(它们应该从红色点到蓝色点).你知道我怎么能改变它吗？如果可能的话,你能告诉我增加箭头大小的方法吗？我尝试了另一个主题(计算和绘图矢量字段)中的建议:

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's 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()

UPDATE

您为要计算的函数添加了规范.它包含依赖于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是系数,并且我计算出CkD1 = sin(2pik)/(k ^ 2 + 1 ^ 2)(这里我使用k和l作为来自n和m之和的索引之一).

我有几个问题是:这两个sin(2*pi*k)和sin(2*pi*k/2)(在两个相互竞争的版本中前因子总是零为整k,给你一个恒定的零V每(x,y)此外,在你的代码,你必须在三角函数,这是从缺少神奇频率的因素.如果你乘以你x的4e-3,你会大幅改变你的功能的空间依赖性(通过将波长改变大约一千倍).所以你应该真正决定你的功能是什么.

所以这是一个解决方案,我假设

V(x,y)= sum_ {k,l = 1到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 + 1 ^ 2))*1E-6

这是您的各种版本的函数的组合,以及sin(2*pi*k/4)在prefactor中的修改,以便具有非零函数.在您找到合适的数学模型后,我希望您能够根据实际需要确定数值因子.

所以这是完整的代码:

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,无论如何它们应该是无限的.

渐变采用与以前相同的方式,但是当转换为numpy函数时,lambdify您必须设置一个额外的字符串参数,'numpy'.这将使得生成的numpy lambda接受数组输入(实际上它将使用np.sin而不是math.sin和相同cos).

我还将网格的定义从更改array为np.linspace:这通常更方便.由于你的函数在整数网格点几乎不变,我创建了一个更密集的网格用于绘图(51点,同时保持原始极限(-10,10)固定).

为清楚起见,我还包括了一些图:a contourf表示函数的值(轮廓线应始终与梯度向量正交),以及一个颜色条以指示函数的值.这是结果:

组成显然不是最好的,但我不想偏离你的规格太多.这个图中的箭头实际上几乎看不到,但正如你所看到的(并且从定义中也很明显V),你的函数是周期性的,所以如果用较小的限制和较少的网格点绘制相同的东西,你会看到更多的特征和更大的箭头.