2015 I ( TA)

Size: px

Start display at page:

Download "2015 I ( TA)"

いおりおとべ
5 years ago
Views:

1 2015 I ( TA)

2 Schrödinger PDE Python u(t, x) x t 2 u(x, t) = k u(t, x) t x2 k k = i h 2m Schrödinger h m 1 ψ(x, t) i h ( 1 ψ(x, t) = i h ) 2 ψ(x, t) t 2m x Cauchy x : x Fourier x x Fourier 2 u(x, t) = k u(t, x) t x2 Fourier F (ξ) = F[F (x)](ξ) 1 2π F (x) = F 1 [F (ξ)](x) 1 2π F (x)e iξx dx F (ξ)e iξx dξ u(x, t) = 1 4kπt e (x x ) 2 4kt u(x, 0)dx 2

3 x 1 4kπt e (x x ) 2 4kt G(x x, t) u(x, t) = G(x x, t)u(x, 0)dx G Green Schrödinger i h ( 1 ψ(x, t) = i h ) 2 ψ(x, t) t 2m x F (p) = F[F (x)](p) 1 2π h F (x) = F 1 [F (p)](x) 1 2π h px i F (x)e h dx F (p)e i px h dp h = 1 ( ) ψ(x, t) = 1 2π h 1 m m(x x ) 2 2π h it ei 2 ht ψ(x, t) = m m(x x ) 2 it ei 2 ht ψ(x, 0)dx = U(x x, t) U(x x, t)ψ(x, 0)dx U(x x, t) G(x x, t) x : x u t=0 ψ t=0 3

4 PDE FTCS ( ) t = 0 t = 0.01 t = 0.02 Python NumPy matplotlib PC equation : 2 tu(x, t) = k x u(t, x) (x : 0 < x < 1, t : 0 < t 0.02, k = 0.5) 2 B.C : u(0, t) = u(1, t) = 1 I.C : u(x, 0) = sin(5πx) + sin(30πx) + sin(80πx) + 1 ======================================================= #-*-coding:utf-8-*- import numpy as np from numpy import linalg as LA import matplotlib.pyplot as plt k=0.5 N=100 dt= step=200 x_ini=0.0 x_end=1.0 dx=(x_end - x_ini)/n dx2=dx**2 r=k*(dt/dx2) def IC(): global u u=np.zeros(n+1) u[0]=1.0 u[n]=1.0 #IC: # for i in range(1,n): # x = x_ini + i*dx 4

5 # if 0<x and x<=0.5: # u[i]=x+1 # elif 0.5<x and x<1: # u[i]=-x+2 #IC: x=np.linspace(x_ini,x_end,n+1) y=np.sin(5*np.pi*x) + np.sin(30*np.pi*x) + np.sin(80*np.pi*x) +1 for i in range(1,n): u[i]=y[i] M=np.linspace(x_ini,x_end,N+1) plt.plot(m,u) def evolve_mat1(): #FTCS A=np.zeros((len(u),len(u)),dtype=float) A[0,0]=A[len(u)-1,len(u)-1]=1 for i in range(1,len(u)-1): A[i,i-1]=r A[i,i]= 1.0-2*r A[i,i+1]=r return A def evolve_mat2(): #Crank-Nicolson myu=r/2 A=np.zeros((len(u),len(u)),dtype=float) B=np.zeros((len(u),len(u)),dtype=float) A[0,0]=A[len(u)-1,len(u)-1]=1 B[0,0]=B[len(u)-1,len(u)-1]=1 for i in range(1,len(u)-1): A[i,i-1]=1 B[i,i-1]=-1 A[i,i]= -2-1/myu B[i,i]= 2-1/myu A[i,i+1]=1 B[i,i+1]=-1 C=np.dot(LA.inv(A),B) return C def evolve(m,b): for i in range(1,step+1): t= i*dt u1=np.dot(m,b) if (i%100)==0 : 5

6 A=np.linspace(x_ini,x_end,N+1) plt.plot(a,u1) b=u1[:] else: plt.legend(["t=0.00","t=0.01","t=0.02"],loc=2) plt.xlabel("x") plt.ylabel("u") plt.axhline(y=0, ls= -, lw=1, color= black ) plt.title("1-d diffusion-eq") plt.show() def main(): IC() Mat=evolve_mat1() # Mat=evolve_mat2() evolve(mat,u) if name == " main ": main() ======================================================= 1: 6

7 FTCS ForwardTimeCenterSpace x i t n u(x, t) = u n i t un i = ( i u n ) i / t 2 x 2 un i = ( u n i+1 2u n i + u n ) i+1 / x 2 i = u n i + k t ( u n x 2 i+1 2u n i + u n ) i 1 1 N N+1 i i : 0 i N i = 0 i = N u = C( ) n u n 0,N i ( k t r ) x i.. N 1 N = r 1 2r r r 1 2r r r 1 2r r 0 0 r 1 2r r 0 1 u n 0 u n 1. u n i.. u n N 1 u n N A = A u n A u u 0,N = A u n + b, u b (i : 1 i N 1) N 1 7

8 (time) (forward) (space) (center) FTCS 1 u 1 u t u t u def evolve_mat1():return A A A=np.zeros((len(u),len(u)),dtype=float) u 0 A A u n 1 Crank-Nicolson 2 3 I FTCS A ( ) 1 u n+1 > u n FTCS r 1 2 Crank-Nicolson VonNeumann 8

9 A max(eigv alue), min(eigv alue) 2: FTCS A 3: Crank-Nicolson A FTCS Crank-Nicolson A 9

10 10

1 matplotlib matplotlib Python matplotlib numpy matplotlib Installing A 2 pyplot matplotlib 1 matplotlib.pyplot matplotlib.pyplot plt import import nu

1 matplotlib matplotlib Python matplotlib numpy matplotlib Installing A 2 pyplot matplotlib 1 matplotlib.pyplot matplotlib.pyplot plt import import nu Python Matplotlib 2016 ver.0.06 matplotlib python 2 3 (ffmpeg ) Excel matplotlib matplotlib doc PDF 2,800 python matplotlib matplotlib matplotlib Gallery Matplotlib Examples 1 matplotlib 2 2 pyplot 2 2.1