# Results blow up when number of intervals is increases (Yee algorithm FDTD, dielectric sphere)

I have been trying to write a program that analyses EM wave scattering by a dielectric sphere for a project.

The reference is Sadiku's book Numerical Methods in electromagnetics Edition 3.

Now the book has code in Matlab but I have to use Python.

I rewrote the code in Python and ran it, but the values Ey/Einc and Ez/Einc blow up, much more so when I increase number of intervals. Even at n = 50 they do not match the results.

Here is my code

'''

import numpy as np
from matplotlib import pyplot as py
import warnings

'''
The followings constants are defined globally, as they are used in the code.
All are in SI units.
'''

pi = np.pi
c = 3e8
e0 = 8.85e-12
mu0 = 4*pi*1e-7

'''
DOCSTRING for Taflove(er, mu, sigma, nu, num)
This funciton is used to calculate some parameters and Taflove's constants.
The inputs are relative permittivity and permeability, conductivity,
frequency of incoming wave and grid number.
Three tuples are returned.
First tuple contains lambd, delta and delta_t parameters.
Second tuple contains R constants.
Third tuple contains C constants.
'''

def Taflove(er, mur, sigma, nu, num):

lambd = c/nu
delta = lambd/num
delta_t = delta/(2*c)

R = delta_t/e0
R_a = (c*delta_t/delta)**2
R_b = delta_t/(mu0*delta)

C_a = (1 - R*(sigma/er))
C_b = R_a/er

return np.array((lambd, delta, delta_t)), np.array((R, R_a, R_b)), np.array((C_a, C_b))

'''
DOCSTRING for
'''

er = (1, 4)          # Relative Permittivity
mur = (1, 1)         # Relative Permeability (1 for non-magnetic)
sigma = (0.1, 0)       # Conductivity (0 for perfect dielectric/lossless medium)
# Conductivity at boundary is taken to be non zero to prevent abrupt results
nu = 2.5e9      # Frequency of EM wave = 2.5 GHz (for Microwave)
num = 40        # Due to computing power

P_list, R_list, C_list1 = Taflove(er, mur, sigma, nu, num)
print(P_list)
print(R_list)
print(C_list1)
input('check1')
_, _, C_list2 = Taflove(er, mur, sigma, nu, num)
C = [C_list1, C_list2]
R, R_a, R_b = R_list
lambd, delta, delta_t = P_list
#1 = air
#2 = sphere
#E = np.array([0, 0, 1])                 # Incoming EM wave is linearly polarized in z direction
#B = np.array([1, 0, 0])/c               # B = k cap cross E/c (in x direction)
#k = 2*pi*np.array([0, 1, 0])/P_list  # k = 2pi/lambda in direction of wave (in y direction)
Size = (19, 39, 19)                     # Grid/Lattice size. More in y because wave is propagating in y driection
Center = (19.5, 20, 19)                 # Center of sphere
WMax = 2                                # Steps in one program looop
TMax =20*5                            # Number of Timesteps
Media = 2                               # Number of Media
JPW = 3                                 # y coordinate of wave

I, J, K = Size
Ex = np.zeros((Size+2, Size+2,Size+2, WMax+1))
Hx = np.zeros((Size+2, Size+2,Size+2, WMax+1))
Ey = np.zeros((Size+2, Size+2,Size+2, WMax+1))
Hy = np.zeros((Size+2, Size+2,Size+2, WMax+1))
Ez = np.zeros((Size+2, Size+2,Size+2, WMax+1))
Hz = np.zeros((Size+2, Size+2,Size+2, WMax+1))
x_ = np.linspace(0, I+1, I+2)
y_ = np.linspace(0, J+1, J+2)
z_ = np.linspace(0, K+1, K+2)
print(x_)
print(y_)
print(z_)
input('check')

pos_func = lambda R : np.sqrt(sum(np.array(R) - np.array(Center))**2)

XMed = np.zeros((Size+2, Size+2,Size+2))
YMed = np.zeros((Size+2, Size+2,Size+2))
ZMed = np.zeros((Size+2, Size+2,Size+2))
for i in range(Size+2):
for j in range(Size + 2):
for k in range(Size+2):
if pos_func((i + 0.5, j, k))<radius_unit:
XMed[i, j, k] = 2
else:
XMed[i, j, k] = 1

YMed[i, j, k] = 2
else:
YMed[i, j, k] = 1

ZMed[i, j, k] = 2
else:
ZMed[i, j, k] = 1

print(XMed)
Ey1 = np.zeros(J+2)
Ez1 = np.zeros(J+2)
print(len(Ey1))

print("Setting up completed...")
input('a')
###########################
t3 = 2
t2 = 1
t1 = 0
for n in range(TMax):
if n%10 == 0:
print(str(100*n/TMax) + " % complete.")

#Applying soft truncation conditons
for k in range(K+1):
for j in range(J+1):
for i in range(I+1):
if i == 0:
if k != K and k!=0:
Hy[0, j, k, t3] = (Hy[1, j, k-1, t1] + Hy[1, j, k, t1] + Hy[1, j, k+1, t1])/3
Hz[0, j, k, t3] = (Hz[1, j, k-1, t1] + Hz[1, j, k, t1] + Hz[1, j, k+1, t1])/3
elif k == K:
Hy[0, j, k, t3] = (Hy[1, j, k-1, t1] + Hy[1, j, k, t1])/2
Hz[0, j, k, t3] = (Hz[1, j, k-1, t1] + Hz[1, j, k, t1])/2
else:
Hy[0, j, k, t3] =(Hy[1, j, k, t1] + Hy[1, j , k+1, t1])/2

Hz[0, j, k, t3] = (Hz[1, j, k, t1] + Hz[1, j, k+1, t1])/2
if j == 0:
Ex[i, 0, k, t3] = Ex[i, 1, k, t1]
Ez[i, 0, k, t3] = Ez[i, 1, k, t1]
elif j == J:
Ex[i, j, k, t3] = Ex[i, j-1, k, t1]
Ez[i, j, k, t3] = Ez[i, j-1, k, t1]
if k == 0:
if i != I and i!= 0:
Ex[i, j, 0, t3] = (Ex[i-1, j, 1, t1] + Ex[i, j, 1, t1] + Ex[i+1, j, 1, t1])/3
Ey[i, j, 0, t3] = (Ey[i-1, j, 1, t2] + Ey[i, j, 1, t2] + Ey[i+1, j, 1, t1])/3
elif i == 0:
Ex[0, j, 0, t3] = (Ex[i, j, 1, t1] + Ex[i+1, j, 1, t1])/2
Ey[0, j, 0, t3] = (Ey[i, j, 1, t1] + Ey[i+1, j, 1, t1])/2
else:
Ex[i, j, 0, t3] = (Ex[i, j, 1, t1] + Ex[i-1, j, 1, t1])/2
Ey[i, j, 0, t3] = (Ey[i, j, 1, t1] + Ey[i-1, j, 1, t1])/2
#Applying Yee's algorithm
Hx[i, j, k, t3] = Hx[i, j, k, t2] + R_b*(Ey[i, j, k+1, t2] - Ey[i, j, k, t2]\
+ Ez[i, j, k, t2] - Ez[i, j+1, k, t2])

Hy[i, j, k, t3] = Hy[i, j, k, t2] + R_b*(Ex[i, j, k, t2] - Ex[i, j, k+1, t2]\
+ Ez[i+1, j, k, t2] - Ez[i, j, k, t2])

Hz[i, j, k, t3] = Hz[i, j, k, t2] + R_b*(Ex[i, j+1, k, t2] - Ex[i, j, k, t2]\
+ Ey[i, j, k, t2] - Ey[i+1, j, k, t2])
if k == K:
Hx[i, j, k, t3] = Hx[i, j, k-1, t3]
Hy[i, j, k, t3] = Hy[i, j, k-1, t3]
if j !=0 and j !=J and k !=0:
M = int(XMed[i, j, k] - 1)
Ex[i, j, k, t3] = C[M]*Ex[i, j, k, t2] + (C[M]/R_list) *(Hz[i, j, k, t3]\
-Hz[i, j-1, k, t3] + Hy[i, j, k-1, t3] - Hy[i, j, k, t3])

if k!=0:
M = int(YMed[i, j, k] - 1)
if i!=0:
Ey[i, j, k, t3] = C[M]*Ey[i, j, k, t2] + (C[M]/R_list) *(Hz[i-1, j, k, t3]\
-Hz[i, j, k, t3] + Hx[i, j, k, t3] - Hx[i, j, k-1, t3])
if Ey[i, j, k, t3]!=0:
#print(Ey[i, j, k, t3])
#input(str(i) +','+ str(j)+',' + str(k))
pass
else:
Ey[i, j, k, t3] = C[M]*Ey[i, j, k, t2] + (C[M]/R_list) *(0\
-Hz[i, j, k, t3] + Hx[i, j, k, t3] - Hx[i, j, k-1, t3])
if j !=0 and j!=J:
M = int(ZMed [i, j, k])-1
if M == 0:
Cam = 1
else:
Cam = C[M]
if i !=0:
Ez [i, j, k,t3] = Cam*Ez[i, j, k, t2]\
+ (C[M]/R_list) *(Hy[i, j, k, t3]\
-Hy[i-1, j, k, t3] + Hx[i, j-1, k, t3] - Hx[i, j, k, t3])

else:
Ez [i, j, k,t3] = Cam*Ez[i, j, k, t2]\
+ (C[M]/R_list) \
*(Hy[i, j, k, t3]\
-0 + Hx[i, j-1, k, t3]- Hx[i, j, k, t3])

if j == JPW:
Ez[i,j,k, t3] = Ez[i,j,k,t3] + np.sin(2*pi*nu*delta_t*(n+1))
if k==K:
Ex[i, j, k + 1, t3] =  Ex[i, j, k , t3]
Ey[i, j, k + 1, t3] =  Ey[i, j, k , t3]

if k == K and n>TMax-num-1:
temp = abs(Ey[I, j, K-1, t3])
if temp>Ey1[j]:
Ey1[j] = temp#/(Ey1[j]+0.0000001)
temp = abs(Ez[I, j, K, t3])
if temp > Ez1[j]:
Ez1[j] =temp#/(Ez1[j]+1e-10)
t1 = int(t2 + 0.)
t2 = int(t3 + 0.)
t3 = int(t3%3)

Integr = np.linspace(1, 100, 100)
print(Ey1)
for i in range(len(Ey1)):
if Ey1[i]<1:
py.plot(Integr[i], Ey1[i], 'r*')
py.xlabel('j')
py.ylabel('Computer Ey/Einc ')
py.grid()
py.show()

print(Ez1)
for i in range(len(Ez1)):
if Ez1[i]<1:
py.plot(Integr[i], Ez1[i], 'r*')
py.xlabel('j')
py.ylabel('Computer Ez/Einc ')
py.grid()
py.show()



Does anyone know what could be the bug? I have been checking this only for the last week. Thanks.

• Potentially troublesome: np.sin(2*pi*nu*delta_t*(n+1)). How big in magnitude can the input to the sine function get there? Generally speaking, for maximum accuracy, you would want sinpi(2*nu ...) instead of sin(2*pi*nu...). Check whether there is sinpi() in Python or a suitable Python library. Otherwise, I would suggest filing a feature request with the Python developers to add sinpi(). Mar 8 at 5:38
• The argument cannot be that big, nu*delta t is of order of ~0.1, n goes to ~100, so the argument is reasonable I think. There is no sinpi() function :( Mar 8 at 5:45
• Probably harmless then. I was thinking this could be the kind of case encountered in this question. I looked over the code and couldn't find anything obviously out of place like an incorrect variable name, or incorrect indexing. Is this a straight port of the MATLAB code in the book, or a heavily modified version? Mar 8 at 6:08
• This is 90% the same version. Should I add the MATLAB code? (it is image not text :( ) Mar 8 at 6:12
• I do not think that is a good idea as it would seem to involve the unauthorized propagation of copyrighted material. Mar 8 at 6:21