薛定谔ampampamp薛定谔方程的pythonmatplotlib实现

　　1 说明：
　　=====
　　1.1 代码来源：感谢原作者，对代码进行修改，注释，删减，bug解决，调试，改进。https://github.com/jakevdp/pySchrodinger  #github地址，别忘记给他点赞
　　1.2 效果图：
　　2 环境：
　　python3.8
　　3 代码分解：
　　========
　　3.1 第1步：导入模块import numpy as np from scipy import fftpack import matplotlib.pyplot as plt from matplotlib import animation
　　3.2 第2步：薛定谔方程，我把它嵌套进入，源代码作为外部库导入class Schrodinger(object):     #初始化     def __init__(self, x, psi_x0, V_x, k0=None, hbar=1, m=1, t0=0.0):         # Validation of array inputs         self.x, psi_x0, self.V_x = map(np.asarray, (x, psi_x0, V_x))         N = self.x.size         assert self.x.shape == (N,)         assert psi_x0.shape == (N,)         assert self.V_x.shape == (N,)          # Validate and set internal parameters         assert hbar > 0         assert m > 0         self.hbar = hbar         self.m = m         self.t = t0         self.dt_ = None         self.N = len(x)         self.dx = self.x[1] - self.x[0]         self.dk = 2 * np.pi / (self.N * self.dx)          # Set momentum scale         if k0 == None:             self.k0 = -0.5 * self.N * self.dk         else:             assert k0 < 0             self.k0 = k0         self.k = self.k0 + self.dk * np.arange(self.N)          self.psi_x = psi_x0         self.compute_k_from_x()          # Variables which hold steps in evolution         self.x_evolve_half = None         self.x_evolve = None         self.k_evolve = None      def _set_psi_x(self, psi_x):         assert psi_x.shape == self.x.shape         self.psi_mod_x = (psi_x * np.exp(-1j * self.k[0] * self.x)                           * self.dx / np.sqrt(2 * np.pi))         self.psi_mod_x /= self.norm         self.compute_k_from_x()      def _get_psi_x(self):         return (self.psi_mod_x * np.exp(1j * self.k[0] * self.x)                 * np.sqrt(2 * np.pi) / self.dx)      def _set_psi_k(self, psi_k):         assert psi_k.shape == self.x.shape         self.psi_mod_k = psi_k * np.exp(1j * self.x[0] * self.dk                                         * np.arange(self.N))         self.compute_x_from_k()         self.compute_k_from_x()      def _get_psi_k(self):         return self.psi_mod_k * np.exp(-1j * self.x[0] * self.dk                                         * np.arange(self.N))      def _get_dt(self):         return self.dt_      def _set_dt(self, dt):         assert dt != 0         if dt != self.dt_:             self.dt_ = dt             self.x_evolve_half = np.exp(-0.5 * 1j * self.V_x                                          / self.hbar * self.dt)             self.x_evolve = self.x_evolve_half * self.x_evolve_half             self.k_evolve = np.exp(-0.5 * 1j * self.hbar / self.m                                     * (self.k * self.k) * self.dt)      def _get_norm(self):         return self.wf_norm(self.psi_mod_x)      psi_x = property(_get_psi_x, _set_psi_x)     psi_k = property(_get_psi_k, _set_psi_k)     norm = property(_get_norm)     dt = property(_get_dt, _set_dt)      def compute_k_from_x(self):         self.psi_mod_k = fftpack.fft(self.psi_mod_x)      def compute_x_from_k(self):         self.psi_mod_x = fftpack.ifft(self.psi_mod_k)      def wf_norm(self, wave_fn):         assert wave_fn.shape == self.x.shape         return np.sqrt((abs(wave_fn) ** 2).sum() * 2 * np.pi / self.dx)      def solve(self, dt, Nsteps=1, eps=1e-3, max_iter=1000):           eps = abs(eps)         assert eps > 0         t0 = self.t         old_psi = self.psi_x         d_psi = 2 * eps         num_iter = 0         while (d_psi > eps) and (num_iter <= max_iter):             num_iter += 1             self.time_step(-1j * dt, Nsteps)             d_psi = self.wf_norm(self.psi_x - old_psi)             old_psi = 1. * self.psi_x         self.t = t0      def time_step(self, dt, Nsteps=1):          assert Nsteps >= 0         self.dt = dt         if Nsteps > 0:             self.psi_mod_x *= self.x_evolve_half             #for num_iter in xrange(Nsteps - 1):  #报错，python2的格式             for num_iter in range(Nsteps - 1):  #python3新版格式                 self.compute_k_from_x()                 self.psi_mod_k *= self.k_evolve                 self.compute_x_from_k()                 self.psi_mod_x *= self.x_evolve             self.compute_k_from_x()             self.psi_mod_k *= self.k_evolve             self.compute_x_from_k()             self.psi_mod_x *= self.x_evolve_half             self.compute_k_from_x()             self.psi_mod_x /= self.norm             self.compute_k_from_x()             self.t += dt * Nsteps
　　3.3 第3步：相关函数定义#高斯波x def gauss_x(x, a, x0, k0):      return ((a * np.sqrt(np.pi)) ** (-0.5)             * np.exp(-0.5 * ((x - x0) * 1. / a) ** 2 + 1j * x * k0))  #傅里叶变换高斯波 def gauss_k(k, a, x0, k0):     ＂＂＂     analytical fourier transform of gauss_x(x), above     ＂＂＂     return ((a / np.sqrt(np.pi)) ** 0.5             * np.exp(-0.5 * (a * (k - k0)) ** 2 - 1j * (k - k0) * x0))  #角度 def theta(x):     x = np.asarray(x)     y = np.zeros(x.shape)     y[x > 0] = 1.0     return y  #障碍物 def square_barrier(x, width, height):     return height * (theta(x) - theta(x - width))
　　3.4 第4步：相关参数定义dt = 0.01 N_steps = 50 t_max = 120 frames = int(t_max / float(N_steps * dt))  # specify constants hbar = 1.0   # planck＂s constant m = 1.9      # particle mass  # specify range in x coordinate N = 2 ** 11 dx = 0.1 x = dx * (np.arange(N) - 0.5 * N)  # specify potential V0 = 1.5 L = hbar / np.sqrt(2 * m * V0) a = 3 * L x0 = -60 * L V_x = square_barrier(x, a, V0) V_x[x < -98] = 1E6 V_x[x > 98] = 1E6  # specify initial momentum and quantities derived from it p0 = np.sqrt(2 * m * 0.2 * V0) dp2 = p0 * p0 * 1. / 80 d = hbar / np.sqrt(2 * dp2)  k0 = p0 / hbar v0 = p0 / m psi_x0 = gauss_x(x, d, x0, k0)  # define the Schrodinger object which performs the calculations S = Schrodinger(x=x,                 psi_x0=psi_x0,                 V_x=V_x,                 hbar=hbar,                 m=m,                 k0=-28)
　　3.5 第5步：matplotlib相关设置#图片初始化定义和参数定义 fig = plt.figure() #图1的x坐标 xlim = (-100, 100) #图2的k坐标==类似x坐标 klim = (-5, 5)  # top axes show the x-space data ymin = 0 ymax = V0 ax1 = fig.add_subplot(211, xlim=xlim,                       ylim=(ymin - 0.2 * (ymax - ymin),                             ymax + 0.2 * (ymax - ymin))) psi_x_line, = ax1.plot([], [], c=＂r＂, label=r＂$|psi(x)|#39;) V_x_line, = ax1.plot([], [], c=＂k＂, label=r＂$V(x)#39;) center_line = ax1.axvline(0, c=＂k＂, ls=＂:＂, label=r＂$x_0 + v_0t#34;)  title = ax1.set_title(＂＂) ax1.legend(prop=dict(size=12)) ax1.set_xlabel(＂$x#39;) ax1.set_ylabel(r＂$|psi(x)|#39;)  # bottom axes show the k-space data ymin = abs(S.psi_k).min() ymax = abs(S.psi_k).max() ax2 = fig.add_subplot(212, xlim=klim,                       ylim=(ymin - 0.2 * (ymax - ymin),                             ymax + 0.2 * (ymax - ymin))) psi_k_line, = ax2.plot([], [], c=＂r＂, label=r＂$|psi(k)|#39;)  p0_line1 = ax2.axvline(-p0 / hbar, c=＂k＂, ls=＂:＂, label=r＂$pm p_0#39;) p0_line2 = ax2.axvline(p0 / hbar, c=＂k＂, ls=＂:＂) mV_line = ax2.axvline(np.sqrt(2 * V0) / hbar, c=＂k＂, ls=＂--＂,                       label=r＂$sqrt{2mV_0}#39;) ax2.legend(prop=dict(size=12)) ax2.set_xlabel(＂$k#39;) ax2.set_ylabel(r＂$|psi(k)|#39;)  V_x_line.set_data(S.x, S.V_x)
　　3.6 第6步：函数初始化和动画设置#第6步：函数初始化和动画设置 def init():     psi_x_line.set_data([], [])     V_x_line.set_data([], [])     center_line.set_data([], [])      psi_k_line.set_data([], [])     title.set_text(＂＂)     return (psi_x_line, V_x_line, center_line, psi_k_line, title)  #动画 def animate(i):     S.time_step(dt, N_steps)     psi_x_line.set_data(S.x, 4 * abs(S.psi_x))     V_x_line.set_data(S.x, S.V_x)     center_line.set_data(2 * [x0 + S.t * p0 / m], [0, 1])      psi_k_line.set_data(S.k, abs(S.psi_k))     title.set_text(＂t = %.2f＂ % S.t)     return (psi_x_line, V_x_line, center_line, psi_k_line, title)
　　3.7 第7步：动画挂起和图片展示anim = animation.FuncAnimation(fig, animate, init_func=init,                                frames=frames, interval=30, blit=True)  #图片展示 plt.show()
　　4 薛定谔：
　　=======
　　4.1 图片：
　　薛定谔，来自奥地利钱币头像
　　4.2 埃尔温·薛定谔：Erwin Schrödinger：
　　（1887年8月12日-1961年1月4日)，男，奥地利物理学家，量子力学奠基人之一，发展了分子生物学。
　　4.3 补充一下：量子力学的研究和应用目前我国是世界第一，真棒！！
　　4.4 薛定谔方程：
　　4.4.1 量子力学中的一个基本方程，也是量子力学的一个基本假定，其正确性只能靠实验来检验。
　　4.4.2 是将物质波的概念和波动方程相结合建立的二阶偏微分方程，可描述微观粒子的运动，每个微观系统都有一个相应的薛定谔方程式，通过解方程可得到波函数的具体形式以及对应的能量，从而了解微观系统的性质。