CODEWITHSUNDEEP
pythonnumpyvectormatplotlib
Ohm
Ohm

Reputation: 2442

Python vector field of ODE of three variables

I am trying to plot a vector field of a ODE model with three variables. I would like to average the vectors along the third axis, and present the vector field together with the information of the standard deviation of their values. The ODE system is:

a  = 1.
b1 = 0.1
b2 = 0.11
c1 = 1.5
c2 = 1.6
d  = 0.75
def dudt(a,b1,b2,u,v1,v2):
    return a*u - b1*u*v1 - b2*u*v2
def dv1dt(d,c1,b1,u,v1):
    return -c1*v1 + d*b1*u*v1
def dv2dt(d,c2,b2,u,v2):
    return -c2*v2 + d*b2*u*v2

The function that I am currently using is:

import numpy as np
import matplotlib.pyplot as plt
def plotVF(mS=None, density= 20,color='k'):
    mB1 = np.linspace(0,1.1,int(density))
    mB2 = np.linspace(0,1.1,int(density))
    if mS==None:
        mS  = np.linspace(0,1.1,int(density))
        B1,B2,S = np.meshgrid(mB1,mB2,mS)
        average=True
    else:
        B1,B2 = np.meshgrid(mB1,mB2)
        S = mS
        average=False
    DB1 = dv1dt(d,c1,b1,S,B1)
    DB2 = dv2dt(d,c2,b2,S,B2)
    DS  = dudt(a,b1,b2,S,B1,B2)
    if average:
        print "Averaging"
        DB1std = np.std(DB1,axis=2)
        DB2std = np.std(DB2,axis=2)
        DB1 = np.mean(DB1,axis=2)
        DB2 = np.mean(DB2,axis=2)
        DS  = np.mean(DS,axis=2)
        vecstd = np.hypot(DB1std,DB2std)
        plt.imshow(vecstd)
        plt.colorbar()
        B1,B2 = np.meshgrid(mB1,mB2)
    M = (np.hypot(DB1, DB2, DS))
    M[ M == 0] = 1.
    DB1=DB1/M
    DB2=DB2/M
    DS=DS/M
    print B1.shape,B2.shape,DB1.shape,DB2.shape
    plt.quiver(B1, B2, DB1, DB2, pivot='mid', color=color)
    plt.xlim(0,1.1), plt.ylim(0,1.1)
    plt.grid('on')
    plt.show()

It gives me that the standard deviation along the third axis is zero, which does not make sense. Vector_Field Someone has an idea what am I doing wrong?

Upvotes: 2

Views: 566

Answers (1)

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

Reputation: 35109

Your code is almost perfectly fine. There's only one problem: you're plotting the colormap with a vanilla call to plt.imshow.

As its name suggests, imshow is used for plotting images. As such, it by default doesn't expect coordinate inputs, just a single array containing the pixel data. This implies that a simple call to imshow will have axes limits corresponding to the number of pixels in your image -- in your case the dimensions of your 2d data arrays. If you take a look at the image directly created by imshow, you'll see that the limits go up to x,y=20. Later you set new limits according to your actual underlying mesh, truncating your plot to the first 2 data points.

The solution is to explicitly tell plt.imshow() where you want your plot to reside in coordinate space:

plt.imshow(vecstd, extent=[B1.min(),B1.max(),B2.min(),B2.max()], origin='lower')

The first keyword argument extent gives the x and y limits into which the data should be plotted. Note the important second keyword argument, origin. By default imshow plots things "upside down" in order to not plot actual images upside down. When you're using imshow to plot stuff defined with Cartesian coordinates, you have to tell it that the origin of the coordinate system should not be the upper left corner of the figure (as for images), bur rather the lower left corner (as for regular plots).

Upvotes: 1

Related Questions