Streamplot in Python Error

Question

I was wondering how I could fix the error in the following code

import numpy as np
import matplotlib.pyplot as plt
from sympy.functions.special.polynomials import assoc_legendre
from scipy.misc import factorial, derivative
import sympy as sym

def main():
    t = 36000
    a=637000000
    H=200
    g=9.81
    x = sym.symbols('x')

    for l in range(1, 6):
        ω=np.sqrt(g*H*l*(l+1))/a

        for n in range(l+1):
            nθ, nφ = 128, 256
            θ, φ = np.linspace(0, np.pi, nθ), np.linspace(0, 2*np.pi, nφ)
            legfun_sym = sym.functions.special.polynomials.assoc_legendre(l, n, x)
            legfun_num = sym.lambdify(x,legfun_sym)

            X, Y = np.meshgrid(θ, φ)

            uθ = (g/(a*ω))*Der_Assoc_Legendre(legfun_num, l, n, X)*np.sin(n*Y-ω*t)
            uφ = (g/(a*ω*np.sin(X)))*Assoc_Legendre(l, n, X)*np.cos(n*Y-ω*t)
            #speed = np.sqrt(uθ**2 + uφ**2)

            fig0, ax = plt.subplots()
            strm = ax.streamplot(φ, θ, uφ, uθ, linewidth=2, cmap=plt.cm.autumn)
            fig0.colorbar(strm.lines)
            plt.show()

def Assoc_Legendre(m, n, X):
    L=[]
    for i in X:
        k=[]
        for j in i:
            k.append(assoc_legendre(m, n, np.cos(j)))
        L.append(k)
    return np.array(L)

def Der_Assoc_Legendre(legfun_num, m, n, X):
    L=[]
    for i in X:
        k=[]
        for j in i:
            k.append(derivative(legfun_num, j, dx=1e-7))
        L.append(k)
    return np.array(L)

if __name__=='__main__':
    main()

The error message 'u' and 'v' must be of shape 'Grid(x,y)' comes up with regard to the strm = ax.streamplot(φ, θ, uφ, uθ, linewidth=2, cmap=plt.cm.autumn) line. How should I fix this?

For reference, I am trying to do a streamplot of $u_{\theta}$ and $u_{\phi}$, where $u_{\theta}=\frac{g}{\omega a}\frac{d}{d\theta}\left(P^n_l\left(cos\theta\right)\right)sin\left(n\phi-\omega t\right)$ and $u_{\phi}=\frac{gn}{\omega a sin\theta}P^n_l\left(cos\theta\right)cos\left(n\phi-\omega t\right)$

EDIT:

This is the current code I have:

import numpy as np
import matplotlib.pyplot as plt
from sympy.functions.special.polynomials import assoc_legendre
from scipy.misc import factorial, derivative
import sympy as sym

def main():
    t = 36000
    a=637000000
    H=200
    g=9.81
    x = sym.symbols('x')
    X, Y = np.mgrid[0.01:np.pi-0.01:100j,0:2*np.pi:100j]

    for l in range(1, 6):
        ω=np.sqrt(g*H*l*(l+1))/a

        for n in range(l+1):
            #nθ, nφ = 128, 256
            #θ, φ = np.linspace(0.001, np.pi-0.001, nθ), np.linspace(0, 2*np.pi, nφ)
            legfun_sym = sym.functions.special.polynomials.assoc_legendre(l, n, x)
            legfun_num = sym.lambdify(x, legfun_sym)
            uθ = (g/(a*ω*np.sin(X)))*Der_Assoc_Legendre(legfun_num, l, n, X)*np.sin(n*Y-ω*t)
            uφ = (g/(a*ω))*Assoc_Legendre(l, n, X)*np.cos(n*Y-ω*t)
            #speed = np.sqrt(uθ**2 + uφ**2)
            fig0, ax = plt.subplots()
            strm = ax.streamplot(Y, X, uθ,uφ, linewidth=0.5, cmap=plt.cm.autumn)
            #fig0.colorbar(strm.lines)
            plt.show()
            print("next")

def Assoc_Legendre(m, n, X):
    L=[]
    for i in X:
        k=[]
        for j in i:
            k.append(assoc_legendre(m, n, np.cos(j)))
        L.append(k)
    return np.float64(np.array(L))

def Der_Assoc_Legendre(legfun_num, m, n, X):
    L=[]
    for i in X:
        k=[]
        for j in i:
            k.append(derivative(legfun_num, j, dx=0.001))
        L.append(k)
    return np.float64(np.array(L))

if __name__=='__main__':
    main()

The current issue seems to be with the derivative function in Der_Assoc_Legendre, that brings up the error ValueError: math domain error after plotting the first plot and onto the second.

Answer 1

While python 3 allows you to use Greek characters as/in variable names, I can assure you that most programmers will find your code unreadable, and it will be a nightmare for others to maintain/develop your code riddled with φ and θ.

Secondly, your code quickly throws a RuntimeWarning concerning a division by zero, which you should most certainly trace down and fix safely.

As for your question, the problem is two-fold. The first problem is that the dimensions of your input don't match in the call to streamline:

>>> print(φ.shape, θ.shape, uφ.shape, uθ.shape)
(256,) (128,) (256, 128) (256, 128)

The trick is that a lot of matplotlib plot functions expect a transposition of its 2d array dimensions, closely related to the weird definition of numpy.meshgrid:

>>> i,j = np.meshgrid(range(3),range(4))
>>> print(i.shape)
(4, 3)

Probably due to this reason, the definition of streamplot is as follows:

 Axes.streamplot(ax, *args, **kwargs)

    Draws streamlines of a vector flow.

    x, y : 1d arrays
        an evenly spaced grid.
    u, v : 2d arrays
        x and y-velocities. Number of rows should match length of y, and the number of columns should match x.

Note the last bit about dimensions. All you need to do is swap x/y, or transpose the angles; you need to check which one of these will lead to a more meaningful plot in your application.

Now, if you fix this, the the following happens:

TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Now this is fishy. All input types should be numeric...right? Well, yeah, but they aren't:

>>> print(φ.dtype, θ.dtype, uφ.dtype, uθ.dtype)
float64 float64 object float64

What's that third object-typed array about?

>>> print(uφ[0,0],uθ[0,0])
+inf -0.00055441014491
>>> print(type(uφ[0,0]),type(uθ[0,0]))
<class 'sympy.core.numbers.Float'> <class 'numpy.float64'>

As @Jelmes noted in a comment. the above type of uφ is the direct consequence of its construction using sympy. If one converts these sympy floats to python or numpy floats before constructing the resulting array, the dtype issue should go away. Whether the remaining infinities (the consequence of division by 0 in 1/sin(X)) will be handled gracefully by streamplot, is another question.