I am newbie in python and doing coding for my physics project which requires to generate a matrix with a variable E

Question

I am newbie in python and doing coding for my physics project which requires to generate a matrix with a variable E for which first element of the matrix has to be solved. Please help me. Thanks in advance. Here is the part of code

import numpy as np
import pylab as pl
import math
import cmath
import sympy as sy
from scipy.optimize import fsolve

#Constants(Values at temp 10K)
hbar = 1.055E-34 
m0=9.1095E-31   #free mass of electron
q= 1.602E-19    
v = [0.510,0,0.510] # conduction band offset in eV
m1= 0.043 #effective mass in In_0.53Ga_0.47As
m2 = 0.072 #effective mass in Al_0.48In_0.52As
d = [-math.inf,100,math.inf]       # dimension of structure in nanometers

'''scaling factor to with units of E in eV, mass in terms of free mass of electron, length in terms 
of nanometers '''
s = (2*q*m0*1E-18)/(hbar)**2 
#print('scaling factor is ',s)
E = sy.symbols('E') #Suppose energy of incoming particle is 0.3eV
m = [0.043,0.072,0.043] #effective mass of electrons in layers
for i in range(3):
       print ('Effective mass of e in layer', i ,'is', m[i])
k=[ ]  #Defining an array for wavevectors in different layers
for i in range(3):   
      k.append(sy.sqrt(s*m[i]*(E-v[i])))
      print('Wave vector in layer',i,'is',k[i])     
    
 x = []  
 for i in range(2):
      x.append((k[i+1]*m[i])/(k[i]*m[i+1]))
      # print(x[i])
 #Define Boundary condition matrix for two interfaces. 
      D0 = (1/2)*sy.Matrix([[1+x[0],1-x[0]], [1-x[0], 1+x[0]]], dtype = complex)
      #print(D0)
  #A = sy.matrix2numpy(D0,dtype=complex)
 D1 = (1/2)*sy.Matrix([[1+x[1],1-x[1]], [1-x[1], 1+x[1]]], dtype = complex)
 #print(D1)   
 #a=eye(3,3)
  #print(a)

  #Define Propagation matrix for 2nd layer or quantum well
    #print(d[1])
    #print(k[1])
     P1 = 1*sy.Matrix([[sy.exp(-1j*k[1]*d[1]), 0],[0, sy.exp(1j*k[1]*d[1])]], dtype = complex)
    #print(P1)
    print("abs")
     T= D0*P1*D1
     #print('Transfer Matrix is given by:',T)
     #print('Dimension of tranfer matrix T is' ,T.shape)

      #print(T[0,0]
      # I want to solve T{0,0} = 0 equation for E
      def f(x):
           return T[0,0]

      x0= 0.5 #intial guess
       x = fsolve(f, x0)
        print("E is",x)
        '''
        y=sy.Eq(T[0,0],0)
         z=sy.solve(y,E)
       print('z',z)
       '''

**The main part i guess is the part of the code where i am trying to solve the equation.***Steps I am following:

1. Defining a symbol E by using sympy
2. Generating three matrices which involves sum formulae and with variable E
3. Generating a matrix T my multiplying those 3 matrices,note that elements are complex and involves square roots of negative number.
4. I need to solve first element of this matrix T[0,0]=0,for variable E and find out value of E. I used fsolve for soving T[0,0]=0.*

Answer 1

Just a note for future questions, please leave out unused imports such as numpy and leave out zombie code like # a = eye(3,3). This helps keep the code as clean and short as possible. Also, the sample code would not run because of indentation problems, so when you copy and paste code, make sure it works before you do so. Always try to make your questions as short and modular as possible.

The expression of T[0,0] is too complex to solve analytically by SymPy so numerical approximation is needed. This leaves 2 options:

1. using SciPy's solvers which are advanced but require type casting to float values since SciPy does not deal with SymPy objects in any way.
2. using SymPy's root solvers which are less advanced but are probably simpler to use.

Both of these will only ever produce a single number as output since you can't expect numeric solvers to find every root. If you wanted to find more than one, then I advise that you use a list of points that you want to use as initial values, input each of them into the solvers and keep track of the distinct outputs. This will however never guarantee that you have obtained every root.

Only mix SciPy and SymPy if you are comfortable using both with no problems. SciPy doesn't play at all with SymPy and you should only have list, float, and complex instances when working with SciPy.

import math
import sympy as sy
from scipy.optimize import newton

# Constants(Values at temp 10K)
hbar = 1.055E-34
m0 = 9.1095E-31  # free mass of electron
q = 1.602E-19
v = [0.510, 0, 0.510]  # conduction band offset in eV
m1 = 0.043  # effective mass in In_0.53Ga_0.47As
m2 = 0.072  # effective mass in Al_0.48In_0.52As
d = [-math.inf, 100, math.inf]  # dimension of structure in nanometers

'''scaling factor to with units of E in eV, mass in terms of free mass of electron, length in terms 
of nanometers '''
s = (2 * q * m0 * 1E-18) / hbar ** 2
E = sy.symbols('E')  # Suppose energy of incoming particle is 0.3eV
m = [0.043, 0.072, 0.043]  # effective mass of electrons in layers
for i in range(3):
    print('Effective mass of e in layer', i, 'is', m[i])
k = []  # Defining an array for wavevectors in different layers
for i in range(3):
    k.append(sy.sqrt(s * m[i] * (E - v[i])))
    print('Wave vector in layer', i, 'is', k[i])

x = []
for i in range(2):
    x.append((k[i + 1] * m[i]) / (k[i] * m[i + 1]))
    # Define Boundary condition matrix for two interfaces.
D0 = (1 / 2) * sy.Matrix([[1 + x[0], 1 - x[0]], [1 - x[0], 1 + x[0]]], dtype=complex)
D1 = (1 / 2) * sy.Matrix([[1 + x[1], 1 - x[1]], [1 - x[1], 1 + x[1]]], dtype=complex)

# Define Propagation matrix for 2nd layer or quantum well
P1 = 1 * sy.Matrix([[sy.exp(-1j * k[1] * d[1]), 0], [0, sy.exp(1j * k[1] * d[1])]], dtype=complex)
print("abs")
T = D0 * P1 * D1

# did not converge for 0.5
x0 = 0.75

# method 1:
def f(e):
    # evaluate T[0,0] at e and remove all sympy related things.
    result = complex(T[0, 0].replace(E, e))
    return result

solution1 = newton(f, x0)
print(solution1)

# method 2:
solution2 = sy.nsolve(T[0,0], E, x0)
print(solution2)

This prints:

(0.7533104353644469-0.023775286117722193j)
1.00808496181754 - 0.0444042144405285*I

Note that the first line is a native Python complex instance while the second is an instance of SymPy's complex number. One can convert the second simply with print(complex(solution2)).

Now, you'll notice that they produce different numbers but both are correct. This function seems to have a lot of zeros as can be shown from the Geogebra plot:

The red axis is Re(E), green is Im(E) and blue is |T[0,0]|. Each of those "spikes" are probably zeros.