In [1]:

    

%matplotlib inline
import numpy as np
import scipy as sp
import pandas as pd
from scipy import interpolate
from scipy import stats
import chaospy as cp

import matplotlib.pyplot as plt


    

In [2]:

    

R = 40.
H = 60
x_w = -40#30.
y_w = 60#90.
Rw = 40#20
seed = 3

def f(x,y,seed):
    np.random.seed(seed)
    return - 3.*np.exp( -((x-x_w)**2. + (y-y_w)**2.)/(Rw**2.) ) + \
           - 3.*np.exp( -((x-30)**2. + (y-40)**2.)/(20**2.) ) + \
           - 1.5*np.exp( -((x-np.random.uniform(-60,60))**2. + (y-np.random.uniform(40,70))**2.)/(60.**2.) ) #+ 4.*(y/H)**0.0 
    
def f1(x,y):
    return - 3.*np.exp( -((x-x_w)**2. + (y-y_w)**2.)/(Rw**2.) )

def f3(x,y):
    return - 3.*np.exp( -((x-30)**2. + (y-40)**2.)/(20**2.) )
    
def f2(x,y,seed):
    np.random.seed(seed)
    return - 1.5*np.exp( -((x-np.random.uniform(-60,60))**2. + (y-np.random.uniform(40,70))**2.)/(60.**2.) ) #+ 4.*(y/H)**0.0


    

In [3]:

    

x_mesh,y_mesh = np.meshgrid(np.arange(-50.,100.,0.1),np.arange(0.,120.,0.1))

F = f(x_mesh,y_mesh,seed)


    

In [4]:

    

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
print


    

In [5]:

    

def MC_int(f,seed):
    x_MC = np.random.uniform(-R,R,9e6)
    y_MC = np.random.uniform(-R,R,9e6) + H

    r = np.sqrt(x_MC**2. + (y_MC-H)**2.)
    ind_in = np.where(r <= R)

    u_eq_MC = np.mean(f(x_MC[ind_in],y_MC[ind_in],seed))
    return u_eq_MC, x_MC[ind_in], y_MC[ind_in]


    

In [6]:

    

u_eq_MC, x_MC, y_MC = MC_int(f,seed)


    

In [7]:

    

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
plt.plot(x_MC,y_MC,'.k',alpha=0.05)
x_plot = np.arange(-50.,50.1,0.1)
print len(x_MC)


    

    




7067113






    

In [8]:

    

def QMC_int(f,seed):
    dist_y = cp.Uniform(-R,R)
    dist_x = cp.Uniform(-(R**2.-dist_y**2.)**0.5,(R**2.-dist_y**2.)**0.5)#cp.Uniform(-R,R)#
    dist = cp.J(dist_y,dist_x)

    sample = dist.sample(size=40,rule='M',antithetic=[1])
    y_MC = sample[0,:]+H
    x_MC = sample[1,:]

    r = np.sqrt(x_MC**2. + (y_MC-H)**2.)
    ind_in = np.where(r <= R)
        
    u_eq_QMC = np.mean(f(x_MC[ind_in],y_MC[ind_in],seed))
    
    return u_eq_QMC, x_MC[ind_in], y_MC[ind_in]


    

In [9]:

    

u_eq_QMC, x_MC, y_MC = QMC_int(f,seed)

print u_eq_QMC
print 100*(u_eq_QMC-u_eq_MC)/u_eq_MC, '%'


    

    




-2.60015706384
-1.49032189306 %

In [10]:

    

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
plt.plot(x_MC,y_MC,'ok',alpha=0.75)
x_plot = np.arange(-50.,50.1,0.1)
print len(x_MC)


    

    




40






    

In [11]:

    

def QMC_int2(f,seed):
    dist_r = cp.Triangle(0,R,R)
    dist_th = cp.Uniform(0,360)
    dist = cp.J(dist_r,dist_th)
    dist_U = cp.J(cp.Uniform(),cp.Uniform())
    
    sample = dist.sample(size=40,rule='M')
    r  = sample[0,:]
    th = sample[1,:]+90

    x_MC = r*np.cos(np.deg2rad(th))
    y_MC = r*np.sin(np.deg2rad(th))+H

    ind_in = np.where(r <= R)
    
    u_eq_QMC = np.mean(f(x_MC[ind_in],y_MC[ind_in],seed))
    
    return u_eq_QMC, x_MC[ind_in], y_MC[ind_in]


    

In [12]:

    

u_eq_QMC2, x_MC, y_MC = QMC_int2(f,seed)

print u_eq_QMC2
print 100*(u_eq_QMC2-u_eq_MC)/u_eq_MC, '%'


    

    




-2.63934519397
-0.00563847197049 %

In [13]:

    

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
plt.plot(x_MC,y_MC,'ok',alpha=0.75)
x_plot = np.arange(-50.,50.1,0.1)
print len(x_MC)


    

    




40






    

In [14]:

    

def PCE_int_cart(f,seed):
    dist_x = cp.Uniform(-40,40)
    dist_y = cp.Uniform(-(R**2.-dist_x**2.)**0.5+H,(R**2.-dist_x**2.)**0.5+H)
    dist = cp.J(dist_x,dist_y)
    dist_U = cp.J(cp.Uniform(),cp.Uniform())
    
    M = 4
    P = cp.orth_ttr(M,dist_U)
    nodes_U, weights = cp.generate_quadrature(M, dist_U, rule='G',sparse=True)
    nodes = dist.inv(dist_U.fwd(nodes_U)) 

    x_MC = nodes[0,:]
    y_MC = nodes[1,:]
    
    sample_u_train = f(x_MC,y_MC,seed)
    u_hat = cp.fit_quadrature(P, nodes, weights, sample_u_train)
    u_eq_PCE = cp.E(u_hat,dist_U)
   
    return u_eq_PCE, x_MC, y_MC, weights


    

In [15]:

    

u_eq_PCE_cart, x_MC, y_MC, weights = PCE_int_cart(f,seed)

print u_eq_PCE_cart
print 100*(u_eq_PCE_cart-u_eq_MC)/u_eq_MC, '%'


    

    




-2.72844356419
3.36994675408 %

In [16]:

    

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
plt.plot(x_MC,y_MC,'ok',alpha=0.75)

print len(x_MC)


    

    




53






    

In [17]:

    

def PCE_int(f,seed):
    dist_r = cp.Triangle(0,R,R)
    dist_th = cp.Uniform(0,360)
    dist = cp.J(dist_r,dist_th)
    dist_U = cp.J(cp.Uniform(),cp.Uniform())
    
    M = 9
    P = cp.orth_ttr(M,dist_U)
    nodes_U, weights = cp.generate_quadrature([3,M], dist_U, rule='G')
    nodes = dist.inv(dist_U.fwd(nodes_U)) 

    r = nodes[0,:]
    th = nodes[1,:]-90
    x_MC = r*np.cos(np.deg2rad(th))
    y_MC = r*np.sin(np.deg2rad(th))+H

    sample_u_train = f(x_MC,y_MC,seed)
    u_hat = cp.fit_quadrature(P, nodes, weights, sample_u_train)
    
    u_eq_PCE = cp.E(u_hat,dist_U)
   
    return u_eq_PCE, x_MC, y_MC, nodes, weights


    

In [18]:

    

u_eq_PCE, x_MC, y_MC, nodes, weights = PCE_int(f,seed)

print u_eq_PCE
print 100*(u_eq_PCE-u_eq_MC)/u_eq_MC, '%'


    

    




-2.64018099787
0.0260268356383 %

In [19]:

    

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
plt.plot(x_MC,y_MC,'ok',alpha=0.75)
print len(x_MC)


    

    




40






    

In [20]:

    

for i in range(len(weights)):
    print '[', '{0:0.20f}'.format(nodes[0,i]/R), ', ', '{0:0.20f}'.format(np.mod(np.deg2rad(nodes[1,i]-90.),2.*np.pi)),  ', ', '{0:0.20f}'.format(weights[i]), '], '
print    

for i in range(len(weights)):
    print '      r_node('+str(i+1)+') =', '{0:0.20f}'.format(nodes[0,i]/R)
    print '      th_node('+str(i+1)+') =', '{0:0.20f}'.format(np.mod(np.deg2rad(nodes[1,i]-90.),2.*np.pi))
    print '      weight('+str(i+1)+') =', '{0:0.20f}'.format(weights[i])
    print

#print nodes[0,:]/R
#print nodes[1,:]
#print weights


    

    




[ 0.26349922998554242692 ,  4.79436403870179805864 ,  0.00579798753740115753 ], 
[ 0.26349922998554242692 ,  5.13630491629471475079 ,  0.01299684397858970851 ], 
[ 0.26349922998554242692 ,  5.71955352542765460555 ,  0.01905256317618122044 ], 
[ 0.26349922998554242692 ,  0.20924454049880022999 ,  0.02341643323656225281 ], 
[ 0.26349922998554242692 ,  1.10309379714216659885 ,  0.02569988335562909190 ], 
[ 0.26349922998554242692 ,  2.03849885644762496284 ,  0.02569988335562912660 ], 
[ 0.26349922998554242692 ,  2.93234811309099407950 ,  0.02341643323656214179 ], 
[ 0.26349922998554242692 ,  3.70522443534172518653 ,  0.01905256317618119616 ], 
[ 0.26349922998554242692 ,  4.28847304447466459720 ,  0.01299684397858971198 ], 
[ 0.26349922998554242692 ,  4.63041392206758217753 ,  0.00579798753740114539 ], 
[ 0.57446451431535072718 ,  4.79436403870179805864 ,  0.01086984853977092380 ], 
[ 0.57446451431535072718 ,  5.13630491629471475079 ,  0.02436599330905551281 ], 
[ 0.57446451431535072718 ,  5.71955352542765460555 ,  0.03571902745281423097 ], 
[ 0.57446451431535072718 ,  0.20924454049880022999 ,  0.04390024659093685194 ], 
[ 0.57446451431535072718 ,  1.10309379714216659885 ,  0.04818117282305908744 ], 
[ 0.57446451431535072718 ,  2.03849885644762496284 ,  0.04818117282305915683 ], 
[ 0.57446451431535072718 ,  2.93234811309099407950 ,  0.04390024659093664378 ], 
[ 0.57446451431535072718 ,  3.70522443534172518653 ,  0.03571902745281418240 ], 
[ 0.57446451431535072718 ,  4.28847304447466459720 ,  0.02436599330905552321 ], 
[ 0.57446451431535072718 ,  4.63041392206758217753 ,  0.01086984853977089951 ], 
[ 0.81852948743000586429 ,  4.79436403870179805864 ,  0.01086984853977090992 ], 
[ 0.81852948743000586429 ,  5.13630491629471475079 ,  0.02436599330905548505 ], 
[ 0.81852948743000586429 ,  5.71955352542765460555 ,  0.03571902745281418934 ], 
[ 0.81852948743000586429 ,  0.20924454049880022999 ,  0.04390024659093679643 ], 
[ 0.81852948743000586429 ,  1.10309379714216659885 ,  0.04818117282305903193 ], 
[ 0.81852948743000586429 ,  2.03849885644762496284 ,  0.04818117282305909438 ], 
[ 0.81852948743000586429 ,  2.93234811309099407950 ,  0.04390024659093658826 ], 
[ 0.81852948743000586429 ,  3.70522443534172518653 ,  0.03571902745281413383 ], 
[ 0.81852948743000586429 ,  4.28847304447466459720 ,  0.02436599330905549199 ], 
[ 0.81852948743000586429 ,  4.63041392206758217753 ,  0.01086984853977088737 ], 
[ 0.96465960618086743494 ,  4.79436403870179805864 ,  0.00579798753740116100 ], 
[ 0.96465960618086743494 ,  5.13630491629471475079 ,  0.01299684397858971545 ], 
[ 0.96465960618086743494 ,  5.71955352542765460555 ,  0.01905256317618123432 ], 
[ 0.96465960618086743494 ,  0.20924454049880022999 ,  0.02341643323656226669 ], 
[ 0.96465960618086743494 ,  1.10309379714216659885 ,  0.02569988335562910925 ], 
[ 0.96465960618086743494 ,  2.03849885644762496284 ,  0.02569988335562914394 ], 
[ 0.96465960618086743494 ,  2.93234811309099407950 ,  0.02341643323656215567 ], 
[ 0.96465960618086743494 ,  3.70522443534172518653 ,  0.01905256317618120657 ], 
[ 0.96465960618086743494 ,  4.28847304447466459720 ,  0.01299684397858972065 ], 
[ 0.96465960618086743494 ,  4.63041392206758217753 ,  0.00579798753740114886 ], 

      r_node(1) = 0.26349922998554242692
      th_node(1) = 4.79436403870179805864
      weight(1) = 0.00579798753740115753

      r_node(2) = 0.26349922998554242692
      th_node(2) = 5.13630491629471475079
      weight(2) = 0.01299684397858970851

      r_node(3) = 0.26349922998554242692
      th_node(3) = 5.71955352542765460555
      weight(3) = 0.01905256317618122044

      r_node(4) = 0.26349922998554242692
      th_node(4) = 0.20924454049880022999
      weight(4) = 0.02341643323656225281

      r_node(5) = 0.26349922998554242692
      th_node(5) = 1.10309379714216659885
      weight(5) = 0.02569988335562909190

      r_node(6) = 0.26349922998554242692
      th_node(6) = 2.03849885644762496284
      weight(6) = 0.02569988335562912660

      r_node(7) = 0.26349922998554242692
      th_node(7) = 2.93234811309099407950
      weight(7) = 0.02341643323656214179

      r_node(8) = 0.26349922998554242692
      th_node(8) = 3.70522443534172518653
      weight(8) = 0.01905256317618119616

      r_node(9) = 0.26349922998554242692
      th_node(9) = 4.28847304447466459720
      weight(9) = 0.01299684397858971198

      r_node(10) = 0.26349922998554242692
      th_node(10) = 4.63041392206758217753
      weight(10) = 0.00579798753740114539

      r_node(11) = 0.57446451431535072718
      th_node(11) = 4.79436403870179805864
      weight(11) = 0.01086984853977092380

      r_node(12) = 0.57446451431535072718
      th_node(12) = 5.13630491629471475079
      weight(12) = 0.02436599330905551281

      r_node(13) = 0.57446451431535072718
      th_node(13) = 5.71955352542765460555
      weight(13) = 0.03571902745281423097

      r_node(14) = 0.57446451431535072718
      th_node(14) = 0.20924454049880022999
      weight(14) = 0.04390024659093685194

      r_node(15) = 0.57446451431535072718
      th_node(15) = 1.10309379714216659885
      weight(15) = 0.04818117282305908744

      r_node(16) = 0.57446451431535072718
      th_node(16) = 2.03849885644762496284
      weight(16) = 0.04818117282305915683

      r_node(17) = 0.57446451431535072718
      th_node(17) = 2.93234811309099407950
      weight(17) = 0.04390024659093664378

      r_node(18) = 0.57446451431535072718
      th_node(18) = 3.70522443534172518653
      weight(18) = 0.03571902745281418240

      r_node(19) = 0.57446451431535072718
      th_node(19) = 4.28847304447466459720
      weight(19) = 0.02436599330905552321

      r_node(20) = 0.57446451431535072718
      th_node(20) = 4.63041392206758217753
      weight(20) = 0.01086984853977089951

      r_node(21) = 0.81852948743000586429
      th_node(21) = 4.79436403870179805864
      weight(21) = 0.01086984853977090992

      r_node(22) = 0.81852948743000586429
      th_node(22) = 5.13630491629471475079
      weight(22) = 0.02436599330905548505

      r_node(23) = 0.81852948743000586429
      th_node(23) = 5.71955352542765460555
      weight(23) = 0.03571902745281418934

      r_node(24) = 0.81852948743000586429
      th_node(24) = 0.20924454049880022999
      weight(24) = 0.04390024659093679643

      r_node(25) = 0.81852948743000586429
      th_node(25) = 1.10309379714216659885
      weight(25) = 0.04818117282305903193

      r_node(26) = 0.81852948743000586429
      th_node(26) = 2.03849885644762496284
      weight(26) = 0.04818117282305909438

      r_node(27) = 0.81852948743000586429
      th_node(27) = 2.93234811309099407950
      weight(27) = 0.04390024659093658826

      r_node(28) = 0.81852948743000586429
      th_node(28) = 3.70522443534172518653
      weight(28) = 0.03571902745281413383

      r_node(29) = 0.81852948743000586429
      th_node(29) = 4.28847304447466459720
      weight(29) = 0.02436599330905549199

      r_node(30) = 0.81852948743000586429
      th_node(30) = 4.63041392206758217753
      weight(30) = 0.01086984853977088737

      r_node(31) = 0.96465960618086743494
      th_node(31) = 4.79436403870179805864
      weight(31) = 0.00579798753740116100

      r_node(32) = 0.96465960618086743494
      th_node(32) = 5.13630491629471475079
      weight(32) = 0.01299684397858971545

      r_node(33) = 0.96465960618086743494
      th_node(33) = 5.71955352542765460555
      weight(33) = 0.01905256317618123432

      r_node(34) = 0.96465960618086743494
      th_node(34) = 0.20924454049880022999
      weight(34) = 0.02341643323656226669

      r_node(35) = 0.96465960618086743494
      th_node(35) = 1.10309379714216659885
      weight(35) = 0.02569988335562910925

      r_node(36) = 0.96465960618086743494
      th_node(36) = 2.03849885644762496284
      weight(36) = 0.02569988335562914394

      r_node(37) = 0.96465960618086743494
      th_node(37) = 2.93234811309099407950
      weight(37) = 0.02341643323656215567

      r_node(38) = 0.96465960618086743494
      th_node(38) = 3.70522443534172518653
      weight(38) = 0.01905256317618120657

      r_node(39) = 0.96465960618086743494
      th_node(39) = 4.28847304447466459720
      weight(39) = 0.01299684397858972065

      r_node(40) = 0.96465960618086743494
      th_node(40) = 4.63041392206758217753
      weight(40) = 0.00579798753740114886

In [446]:

    

r_aux = nodes[0,:]/R
th_aux = np.mod(np.deg2rad(nodes[1,:]-90.),2.*np.pi)

fig = plt.figure(figsize=[4,4])
ax = fig.add_subplot(111)
ax.scatter(r_aux*np.cos(th_aux),r_aux*np.sin(th_aux),s=30,c=weights,cmap=plt.cm.gray_r)
circ = plt.Circle([0.,0.],1.,ec='k',fill=False)
ax.add_patch(circ)
ax.set_xlim([-1.1,1.1])
ax.set_ylim([-1.1,1.1])
ax.set_axis_off()
fig.savefig('PCE_Quad_rotor.pdf')


    

In [276]:

    

M = 9
P = cp.orth_ttr(M,dist_U)
nodes_U, weights = cp.generate_quadrature([3,M], dist_U, rule='G')
nodes = dist.inv(dist_U.fwd(nodes_U)) 

r = nodes[0,:]
th = nodes[1,:]+90
x_MC = r*np.cos(np.deg2rad(th))
y_MC = r*np.sin(np.deg2rad(th))+H

sample_u_train_1 = f1(x_MC,y_MC)
u_hat = cp.fit_quadrature(P, nodes, weights, sample_u_train_1)
u_eq_PCE_1 = cp.E(u_hat,dist_U)

sample_u_train_2 = f2(x_MC,y_MC,seed)
u_hat = cp.fit_quadrature(P, nodes, weights, sample_u_train_2)
u_eq_PCE_2 = cp.E(u_hat,dist_U)

sample_u_train_3 = f3(x_MC,y_MC)
u_hat = cp.fit_quadrature(P, nodes, weights, sample_u_train_3)
u_eq_PCE_3 = cp.E(u_hat,dist_U)


    

In [277]:

    

u_eq_PCE = u_eq_PCE_1 + u_eq_PCE_2 + u_eq_PCE_3

print u_eq_PCE
print 100*(u_eq_PCE-u_eq_MC)/u_eq_MC, '%'


    

    




-2.62656050526
-0.489999815812 %

In [290]:

    

def gaussN(R, func, varargin, NGr=4, NGth=4):
    """Calculate numerically the gauss integration.
    [1] eq. 38

    Inputs
    ----------
    R (float): Wind turbine radius [m]
    func (function): Wind speed function
    varargin: Other arguments for the function besides [r,te]
    NG (int): Number of Ga

    Outputs
    ----------
    Ua (float):
    """
    A = np.pi*R**2
    #coefficients
    if  (NGr==4)&(NGth==4): #for speed give the full values
        rt = np.array([[ -0.339981043584856, -0.861136311594053,
            0.339981043584856, 0.861136311594053]])
        te = rt.T
        w = np.array([[0.652145154862546, 0.347854845137454,
            0.652145154862546, 0.347854845137454]])
        wt=w
    else:
        rt,w = np.polynomial.legendre.leggauss(NGr)
        rt = np.array([rt])
        #te = rt.T
        w = np.array([w])
        te,wt = np.polynomial.legendre.leggauss(NGr)
        te = np.array([te])
        wt = np.array([wt])
        
    return np.sum((np.pi/4.0)*(R**2./A)*w*wt.T*func(R*(rt+1.0)/2.0,
        np.pi*(te+1.0),*varargin)*(rt+1.0))


    

In [289]:

    

def f_test(x,y,seed):
    np.random.seed(seed)
    return - np.random.uniform(0.5,3.)*np.exp( -((x-np.random.uniform(-60,60))**2. + (y-np.random.uniform(40,70))**2.)/(np.random.uniform(10,80)**2.) ) + \
           - np.random.uniform(0.5,3.)*np.exp( -((x-np.random.uniform(-60,60))**2. + (y-np.random.uniform(40,70))**2.)/(np.random.uniform(10,80)**2.) ) + \
           - np.random.uniform(0.5,3.)*np.exp( -((x-np.random.uniform(-60,60))**2. + (y-np.random.uniform(40,70))**2.)/(np.random.uniform(10,80)**2.) ) 
    
def plot(f,seed):
    x_mesh,y_mesh = np.meshgrid(np.arange(-50.,100.,0.1),np.arange(0.,120.,0.1))

    F2 = f(x_mesh,y_mesh,seed)

    plt.contourf(x_mesh,y_mesh,F2,20)
    cb = plt.colorbar()
    circ = plt.Circle([0.,H],R,ec='k',fill=False)
    ax = plt.gca()
    ax.add_patch(circ)
    print  
    
f_r = lambda r,th,seed: f_test(-r*np.sin(th),r*np.cos(th)+H,seed)


    

In [315]:

    

rel_error_PCE = []
rel_error_QMC = []
rel_error_QMC2 = []
rel_error_GCL = []

seeds = range(1000)

for seed in seeds:  
    u_eq_MC, x_MC, y_MC = MC_int(f_test,seed)
    u_eq_PCE, x_MC, y_MC, nodes, weights = PCE_int(f_test,seed)
    u_eq_QMC, x_MC, y_MC = QMC_int(f_test,seed)
    u_eq_QMC2, x_MC, y_MC = QMC_int2(f_test,seed)
    u_eq_GCL = gaussN(R,f_r,[seed])

    rel_error_PCE = np.append(rel_error_PCE,100*(u_eq_PCE-u_eq_MC)/u_eq_MC)
    rel_error_QMC = np.append(rel_error_QMC,100*(u_eq_QMC-u_eq_MC)/u_eq_MC)
    rel_error_QMC2 = np.append(rel_error_QMC2,100*(u_eq_QMC2-u_eq_MC)/u_eq_MC)
    rel_error_GCL = np.append(rel_error_GCL,100*(u_eq_GCL-u_eq_MC)/u_eq_MC)


    

In [330]:

    

bins = np.linspace(-6,10,50)

fig=plt.figure(figsize=[5,5])
ax = fig.add_subplot(111)
ax.hist(rel_error_PCE,bins=bins,normed=True,histtype='step',label='PCE',color='k')
ax.hist(rel_error_QMC2,bins=bins,normed=True,histtype='step',label='QMC Hammersley',color='b')
ax.hist(rel_error_QMC,bins=bins,normed=True,histtype='step',label='QMC Antithetic',color='r')
ax.hist(rel_error_GCL,bins=bins,normed=True,histtype='step',label='GCL',color='g')
ax.set_xlabel('Relative error in rotor avg multiple wake deficit [%]')
ax.set_ylabel('PDF [-]')
plt.legend(loc=1)
fig.savefig('Rotor_avg_int_approx_comparison.pdf')
print


    

In [ ]: