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.

Andras Deak -- Слава Україні · Accepted Answer · 2017-07-02 22:09:21Z

2

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.

edited Jul 2, 2017 at 22:09

answered Jul 2, 2017 at 21:12

Andras Deak -- Слава Україні

35.4k13 gold badges94 silver badges118 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Jelmes Over a year ago

Although the division by np.sin(X) at 0 and $\pi$ are not great, it seems to me that the difference in dtype steps from the function Assoc_Legendre and Der_Assoc_Legendre. I am not quite sure how to fix this

Andras Deak -- Слава Україні Over a year ago

@Jelmes you're right, I somehow missed that only the problematic one was calculated with sympy. So instead append(float(assoc_legendre(...))), but then you'll still have the infinities (we'll see how streamplot handles those).

Jelmes Over a year ago

If you see the edited post, it gives a reference what in fact I am trying to do. The latex hasn't printed, which is annoying. Possibly caused by an error? anyway, hopefully you can read it :)

Collectives™ on Stack Overflow

Streamplot in Python Error

1 Answer 1

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related