我是 python 的新手，正在为我的物理项目编写代码，该项目需要生成一个带有变量 E 的矩阵

Question

我是 python 的新手，正在为我的物理项目编写代码，该项目需要生成一个带有变量 E 的矩阵，必须求解该矩阵的第一个元素。请帮我。提前致谢。这是部分代码

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)
       '''

**我猜的主要部分是我试图求解方程的代码部分。***我遵循的步骤：

使用sympy定义符号E
生成三个涉及求和公式且变量为E的矩阵
将这 3 个矩阵相乘生成矩阵 T，请注意元素是复数并且涉及负数的平方根。
我需要求解这个矩阵的第一个元素 T[0,0]=0，对于变量 E 并找出 E 的值。我使用 fsolve 来求解 T[0,0]=0。*

Answer 1

只是为以后的问题做一个注释，请省略 numpy 等未使用的导入并省略 # a = eye(3,3) 等僵尸代码。这有助于保持代码尽可能简洁。此外，由于缩进问题，示例代码不会运行，因此当您复制和粘贴代码时，请确保在执行此操作之前它可以正常工作。始终尽量使您的问题尽可能简短和模块化。

T[0,0] 的表达式太复杂，无法通过 SymPy 解析求解，因此需要数值近似。这留下了 2 个选项：

使用 SciPy 的高级求解器，但需要将类型转换为 float 值，因为 SciPy 不以任何方式处理 SymPy 对象。
使用 SymPy 的根求解器，它不太先进但可能更易于使用。

这两个都只会产生一个数字作为输出，因为您不能指望数字求解器找到每个根。如果您想找到多个，那么我建议您使用要用作初始值的点列表，将它们中的每一个输入到求解器中并跟踪不同的输出。但是，这不能保证您已获得所有根。

只有在您可以毫无问题地使用 SciPy 和 SymPy 的情况下才混合使用它们。 SciPy 根本不能与 SymPy 一起玩，在使用 SciPy.[=25= 时你应该只有 list、float 和 complex 个实例]

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)

这会打印：

(0.7533104353644469-0.023775286117722193j)
1.00808496181754 - 0.0444042144405285*I

请注意，第一行是原生 Python complex 实例，而第二行是 SymPy 复数的实例。可以简单地用 print(complex(solution2)).

转换第二个

现在，您会注意到它们产生的数字不同，但都是正确的。这个函数似乎有很多零，从 Geogebra 图中可以看出：

红轴是Re(E)，绿轴是Im(E)，蓝轴是|T[0,0]|。这些“尖峰”中的每一个都可能是零。

我是 python 的新手，正在为我的物理项目编写代码，该项目需要生成一个带有变量 E 的矩阵

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

python

matrix

sympy

equation-solving