How to plot grad(f(x,y))?
Question: I want to calculate and plot a gradient of any scalar function of two
variables. If you really want a concrete example, lets say f=x^2+y^2 where x
goes from -10 to 10 and same for y. How do I calculate and plot grad(f)? The
solution should be vector and I should see vector lines. I am new to python so
please use simple words.

EDIT:

@Andras Deak: thank you for your post, i tried what you suggested and instead
of your test function (fun=3*x^2-5*y^2) I used function that i defined as
V(x,y); this is how the code looks like but it reports an error

    
    
    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'
    

And lets say I remove sin, I get another error:

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

I read tutorial for sympy and it says "The real power of a symbolic
computation system such as SymPy is the ability to do all sorts of
computations symbolically". I get this, I just dont get why I cannot multiply
x and y symbols with float numbers.

What is the way around this? :( Help please!

**UPDATE**

@Andras Deak: I wanted to make things shorter so I removed many constants from
the original formulas for V(x,y) and Cn*Dm. As you pointed out, that caused
the sin function to always return 0 (i just noticed). Apologies for that. I
will update the post later today when i read your comment in details. Big
thanks!

**UPDATE 2** I changed coefficients in my expression for voltage and this is
the result:

[![plot](http://i.stack.imgur.com/FlDLH.png)](http://i.stack.imgur.com/FlDLH.png)

It looks good except that the arrows point in the opposite direction (they are
supposed to go out of the reddish dot and into the blue one). Do you know how
I could change that? And if possible, could you please tell me the way to
increase the size of the arrows? I tried what was suggested in another topic
([Computing and drawing vector
fields](http://stackoverflow.com/questions/25342072/computing-and-drawing-
vector-fields)):

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

This plots only every third arrow and matplotlib does the autoscale but it
doesnt work for me (nothing happens when i add this, for any number that i
enter) You were already of huge help , i cannot thank you enough!


Answer: Here's a solution using `sympy` and `numpy`. This is the first time I use
sympy, so others will/could probably come up with much better and more elegant
solutions.

    
    
    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)
    

now you can use `numgradfun(1,3)` to compute the gradient at `(x,y)==(1,3)`.
This function can then be used for plotting, which you said you can do.

For plotting, you can use, for instance, `matplotlib`'s `quiver`, like so:

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

[![enter image description
here](http://i.stack.imgur.com/iT1Bv.png)](http://i.stack.imgur.com/iT1Bv.png)

## UPDATE

You added a specification for your function to be computed. It contains the
product of terms depending on `x` and `y`, which seems to break my above
solution. I managed to come up with a new one to suit your needs. However,
your function seems to make little sense. From your edited question:

    
    
    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)
    

On the other hand, from your corresponding comment to this answer:

> 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).

I have several problems with this: both `sin(2*pi*k)` and `sin(2*pi*k/2)` (the
two competing versions in the prefactor are always zero for integer `k`,
giving you a constant zero `V` at every `(x,y)`. Furthermore, in your code you
have magical frequency factors in the trigonometric functions, which are
missing from the comment. If you multiply your `x` by `4e-3`, you
_drastically_ change the spatial dependence of your function (by changing the
wavelength by roughly a factor of a thousand). So you should really decide
what your function is.

So here's a solution, where I assumed

> 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

This is a combination of your various versions of the function, with the
modification of `sin(2*pi*k/4)` in the prefactor in order to have a non-zero
function. I expect you to be able to fix the numerical factors to your actual
needs, after you figure out the proper mathematical model.

So here's the full code:

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

I defined your `V(x,y)` function using some auxiliary functions for
transparence. I left the summation cut-offs as literal parameters, `kmax` and
`lmax`: in your code these were 3, in your comment they were said to be 10,
and anyway they should be infinity.

The gradient is taken the same way as before, but when converting to a numpy
function using `lambdify` you have to set an additional string parameter,
`'numpy'`. This will alow the resulting numpy lambda to accept array input
(essentially it will use `np.sin` instead of `math.sin` and the same for
`cos`).

I also changed the definition of the grid from `array` to `np.linspace`: this
is usually more convenient. Since your function is almost constant at integer
grid points, I created a denser mesh for plotting (51 points while keeping
your original limits of (-10,10) fixed).

For clarity I included a few more plots: a `contourf` to show the value of the
function (contour lines should always be orthogonal to the gradient vectors),
and a colorbar to indicate the value of the function. Here's the result:

[![enter image description
here](http://i.stack.imgur.com/juC87.png)](http://i.stack.imgur.com/juC87.png)

The composition is obviously not the best, but I didn't want to stray too much
from your specifications. The arrows in this figure are actually hardly
visible, but as you can see (and also evident from the definition of `V`) your
function is periodic, so if you plot the same thing with smaller limits and
less grid points, you'll see more features and larger arrows.

