In [16]:
import matplotlib 
import matplotlib.pyplot as plt
import glob
import json
import pickle
import ternary
#import dm_tools
import numpy as np
import seaborn as sns
from scipy import stats
sns.set(style="ticks")
import statsmodels.api as sm
import pandas as pd
from scipy.stats import chisquare
from matplotlib import rcParams
from scipy.interpolate import griddata
import math
import matplotlib.cm as cm
import matplotlib.animation as animation

%matplotlib inline



print('Done')
Done
In [17]:
#Heat map generation

def make_plot(c,eh,ax,fig,npts=100):
  x=c[0,:]
  y=c[1,:]
  X,Y=np.meshgrid(np.linspace(x.min(),x.max(),npts),
                  np.linspace(y.min(),y.max(),npts) ) 
  eH=griddata((x,y),eh,(X,Y),method='linear',fill_value=np.mean(eh))
  CM=ax.imshow(np.flipud(eH),cmap=cm.RdYlGn,
            extent=[x.min(),x.max(),y.min(),y.max()],
            interpolation='bicubic')
  fig.colorbar(CM)
  ax.set_xlabel(r'$x$', fontsize=fs+10)
  ax.set_ylabel(r'$y$', fontsize=fs+10)
  return CM

print('Done')
Done
In [48]:
#Perform relaxation method to solve Poisson's equation in 2D

##############################################################################################

#Produce arrays that described the fixed boundary conditions along each edge of the grid
def lx0_boundary(L, h):
    bound=[]
    for k in range(0,L+1):
        x=k*h
#        value=-x**2 #This is where you type in the function that defines the boundary for the left-hand x-axis
        value=+1
        bound.append(value)
    return bound
    
def lx1_boundary(L, h):
    bound=[]
    for k in range(0,L+1):
        x=k*h
#        value=-(1-x)**2 #This is where you type in the function that defines the boundary for the right-hand x-axis
        value=-1
        bound.append(value)
    return bound

def ly0_boundary(L, h):
    bound=[]
    for k in range(0,L+1):
        x=k*h
#        value=-x**2 #This is where you type in the function that defines the boundary for the lower y-axis
        value=0
        bound.append(value)
    return bound

def ly1_boundary(L, h):
    bound=[]
    for k in range(0,L+1):
        x=k*h
#        value=-(1-x)**2 #This is where you type in the function that defines the boundary for the upper y-axis
        value=0
        bound.append(value)
    return bound

#############################################################################################

#Compute the source function over the integration region
def source_function(region, h):
  L=len(region[0])
  sigma=1
  source_region=np.zeros((L,L))
  for k in range(0, L):
    y=k*h
    for l in range(0, L):
        x=l*h
#        f=-10*np.exp((-(x-0.5)**2-(y-0.5)**2)/(2*sigma**2))   #This line defines the source function
#        f=0   #Setting f=0 means you're solving the Laplace equation
#        f=100*y
#        f=0.1*(np.exp(10*y)-1)
        f=0.1*(np.exp(10*y)+np.exp(10*x))
    
        source_region[k][l]=-1*f
  return source_region
        
############################################################################################
    
#Function to generate the NxN region that is to be relaxed, including the boundary conditions along the edges
def region_generator(l0, l1, l2, l3, h):
  L=len(l0)
  region=np.zeros((L,L))
  for k in range(0,L):
    for l in range(0,L):
        region[k][l]=np.random.uniform(-2.0,2.0)
#        region[k][l]=0
  c00=0.5*(l0[0]+l2[0])     #Lower left-hand corner
  c01=0.5*(l1[0]+l2[L-1])   #Lower right-hand corner
  c10=0.5*(l0[L-1]+l3[0])   #Upper left-hand corner
  c11=0.5*(l1[L-1]+l3[L-1]) #Upper right-hand corner

  if l0[0]==l2[0]:
        print('Conditions continuous in lower left-hand corner.')
  else:
        print('Conditions discontinuous in lower left-hand corner')
        
  if l1[0]==l2[L-1]:
        print('Conditions continuous in lower right-hand corner.')
  else:
        print('Conditions discontinuous in lower right-hand corner')
        
  if l0[L-1]==l3[0]:
        print('Conditions continuous in upper left-hand corner.')
  else:
        print('Conditions discontinuous in upper left-hand corner')
        
  if l1[L-1]==l3[L-1]:
        print('Conditions continuous in upper right-hand corner.')
  else:
        print('Conditions discontinuous in upper right-hand corner')
        
        
  region[0]=l2
  region[L-1]=l3
  for i in range(0,L):
    region[i][0]=l0[i]
    region[i][L-1]=l1[i]
    
  region[0][0]=c00      #Write-in the corners of the grid by hand
  region[0][L-1]=c01
  region[L-1][0]=c10
  region[L-1][L-1]=c11
  return region

#######################################################################################

#Function to perform a single relaxation iteration on the grid
def relax_grid(region, source, h):
    L=len(region)
    
    new_region=np.zeros((L,L))
    
    new_region[0]=region[0]
    new_region[L-1]=region[L-1]
    
    for i in range(0,L):
      new_region[i][0]=region[i][0]
      new_region[i][L-1]=region[i][L-1]

    
    for k in range(1,L-1):
        for l in range(1,L-1):
            new_region[k][l]=0.25*(region[k][l+1]+region[k+1][l]+region[k][l-1]+region[k-1][l]-(h**2)*source[k][l])
    
    return new_region

######################################################################################

#Set-up the grid over which the relaxation will be performed

h=0.05      #Define step size

ylim=1      #Upper y-axis limit for the integration region
xlim=ylim   #Right x-axis limit for the integration region

N = np.int(ylim/h)  #Calculate number of grid points
    
lx0 = lx0_boundary(N, h)    #Left-hand x-axis boundary condition
lx1 = lx1_boundary(N, h)    #Right-hand x-axis boundary condition
ly0 = ly0_boundary(N, h)    #Lower y-axis boundary condition
ly1 = ly1_boundary(N, h)    #Upper y-axis boundary condition

region=region_generator(lx0, lx1, ly0, ly1, h)

N=len(region[0])
fs=15
shift=0.05

step_num=50   #Number of relaxation iterations

source=source_function(region, h)

##################################################################################################

#Perform the loop that does the relaxation iteration steps
for k in range(10, step_num+10):

  grid_holder={'x':[], 'y':[], 'solution':[], 'source':[]}
  for i in range(0,N):
    for j in range(0,N):
      grid_holder['x'].append(j*h)
      grid_holder['y'].append(i*h)
      grid_holder['solution'].append(region[i][j])
      grid_holder['source'].append(source[i][j])
  panda_grid=pd.DataFrame(grid_holder)


  plt.figure(figsize=(3,3))
  fig=plt.figure(figsize=(3,3))
  ax=plt.gca()
  slice=make_plot(np.array([panda_grid['x'].values,
                    panda_grid['y'].values]),
                    panda_grid['solution'],ax,fig)
  ax.set_xlim([0-shift,xlim+shift])
  ax.set_ylim([0-shift,ylim+shift])
#  ax.set_xticklabels([0, 0.25, 0.5, 0.75, 1])
#  ax.set_yticklabels([0, 0.25, 0.5, 0.75, 1])
  ax.set_title('Solution, Iteration Number: '+np.str(k-9), fontsize=fs-5)
  plt.savefig('relaxation_iteration_'+np.str(k)+'.png')

  region=relax_grid(region, source, h)
    
############################################################################################
  

print('Done')
Conditions discontinuous in lower left-hand corner
Conditions discontinuous in lower right-hand corner
Conditions discontinuous in upper left-hand corner
Conditions discontinuous in upper right-hand corner
Done
<matplotlib.figure.Figure at 0x8ec37d0>
<matplotlib.figure.Figure at 0xc5994d0>
<matplotlib.figure.Figure at 0xd216f50>
<matplotlib.figure.Figure at 0x119206d0>
<matplotlib.figure.Figure at 0xfafe050>
<matplotlib.figure.Figure at 0x1acc8210>
<matplotlib.figure.Figure at 0x12ad1d10>
<matplotlib.figure.Figure at 0xea6b690>
<matplotlib.figure.Figure at 0xb3e36d0>
<matplotlib.figure.Figure at 0x131fa790>
<matplotlib.figure.Figure at 0x11a77d50>
<matplotlib.figure.Figure at 0x10b6dc90>
<matplotlib.figure.Figure at 0x1b225090>
<matplotlib.figure.Figure at 0x12ad31d0>
<matplotlib.figure.Figure at 0x762c110>
<matplotlib.figure.Figure at 0x762cd10>
<matplotlib.figure.Figure at 0x95fbc50>
<matplotlib.figure.Figure at 0x125811d0>
<matplotlib.figure.Figure at 0x122d4b90>
<matplotlib.figure.Figure at 0x852d450>
<matplotlib.figure.Figure at 0x8d66d50>
<matplotlib.figure.Figure at 0x11917f90>
<matplotlib.figure.Figure at 0x1a77cdd0>
<matplotlib.figure.Figure at 0x5eb9690>
<matplotlib.figure.Figure at 0x122d7d50>
<matplotlib.figure.Figure at 0xf882b10>
<matplotlib.figure.Figure at 0x119398d0>
<matplotlib.figure.Figure at 0x6715650>
<matplotlib.figure.Figure at 0xf86cc50>
<matplotlib.figure.Figure at 0xfc65490>
<matplotlib.figure.Figure at 0xfc6e210>
<matplotlib.figure.Figure at 0xafd4bd0>
<matplotlib.figure.Figure at 0x11372750>
<matplotlib.figure.Figure at 0x86c58d0>
<matplotlib.figure.Figure at 0x1135c910>
<matplotlib.figure.Figure at 0x1135ac10>
<matplotlib.figure.Figure at 0x11d39c90>
<matplotlib.figure.Figure at 0xd648e50>
<matplotlib.figure.Figure at 0x130724d0>
<matplotlib.figure.Figure at 0xd63e510>
<matplotlib.figure.Figure at 0x16b97a10>
<matplotlib.figure.Figure at 0x1007a690>
<matplotlib.figure.Figure at 0x193e63d0>
<matplotlib.figure.Figure at 0x918a110>
<matplotlib.figure.Figure at 0xafddc90>
<matplotlib.figure.Figure at 0x1007d2d0>
<matplotlib.figure.Figure at 0xe66e3d0>
<matplotlib.figure.Figure at 0xee8d610>
<matplotlib.figure.Figure at 0xcf43690>
<matplotlib.figure.Figure at 0x126f8c50>