gistfile1.txt

import numpy as np
import matplotlib.pyplot as plt
from math import *   
import pylab 
from matplotlib import mlab  

dlin_p=3; glub=-1; lamda=0.8; kv=2*3.14/lamda;

dx=lamda/0; dt=dx; dlin=round(23.0/dx);  X_ros=0; 

time=0; t=0; time_max=3000000; X_center=round(dlin/2.0)*dx;  h=1.0 ; m=1860.0; 


print("kolich shagov = ",dlin)
print("lamda = ",lamda)
print("dx,dt = ",dt)
print("k = ",kv)


def f(x):
    x = i*dx - X_center
    if(x>=-dlin_p and x<=dlin_p):
        f=glub
    else:
            f=0.0
    return f 


def Psi(x):
    x = i*dx - X_center
    Psi=e**(-(x-X_ros)**2+kv*x*1j)                  
    return Psi 

U=[]

for i in range(dlin+1):                                                                                         
    x = i*dx - X_center                                                                                         
    U+=[(Psi(i*dx - X_center))]

r=[]

for i in range(1,dlin+1):                                                                                                                                                                         
    r+=[i*dx - X_center] 

l=0
for n in range(1,dlin-1):
    l+=(abs(U[n])**2 + abs(U[n+1])**2)/2



B=1.0+1j*dt/(2.0*m*dx**2) + 1j*f(x)*dt/2.0
C=A=1j*dt/(4.0*m*dx**2)

time1=199 
g=[]

pylab.ion()
#d=[]


while True: 

    time+=1 
    V=[]
    p=[]; p+=[0.0+0j];
    q=[]; q+=[0.0+0j];
    for i in range(1,dlin):
        x = i*dx - X_center
        D=(U[i] + (U[i+1]-2*U[i]+U[i-1])*1j*dt/(4*m*dx**2)-1j*dt*f(x)*U[i]/2)
        w2=B-C*p[i-1]
        p+=[A/w2]
        q+=[(D+A*q[i-1])/w2]
    V+=[0.0+0j]
    k=0.0+0j
    for i in range(dlin-1,0,-1):
        k=p[i]*k+q[i]
        V+=[k]
    V+=[0.0+0j]
    del U
    U=[]
    for i in range(dlin+1):
        U+=[V[dlin-i]] 
    del p; del q; del V;
    time1+=1
    if (time1==200):
        del g
        #del d
        g=[]
        l=0
        #d=[]
        for n in range(1,dlin-1):
            l+=(abs(U[n])**2 + abs(U[n+1])**2)/2
 
        for i in range(1,dlin+1):
            g+=[abs(U[i])**2]#[U[i].real]
            #d+=[abs(U[i])**2]
        time1=0
        pylab.clf()
        plt.plot([-dx*dlin/2,-dlin_p,-dlin_p,dlin_p,dlin_p,dx*dlin/2],[0,0,glub,glub,0,0])
        #kve=np.array(d)
        kx=np.array(r)
        ky=np.array(g)
        pylab.plot(kx,ky)
        #pylab.plot(kx,kve)
        plt.grid(True)
        plt.text(-9,1.2, 'k = ' );plt.text(-8,1.2, round(kv,2));plt.text(-5,1.2, 'psi = ' );plt.text(-3.5,1.2, round(l*dx,2));plt.text(0,1.2, '  lamda = ' );plt.text(3,1.2, lamda);plt.text(5,1.2, 't =' );plt.text(6.1,1.2, round(time*dt))
        plt.axis([-11,11, -1.1, 1.1])
        pylab.draw()
        #fname = '_tmp%03d.png'%time
        #plt.savefig(fname)

    if(time==time_max):                                        
        break
pylab.close()  










import numpy as np
import matplotlib.pyplot as plt
from math import *   
import pylab 
from matplotlib import mlab  

dlin_p=1.1; glub=-0.5; lamda=0.3; kv=2*3.14/lamda;

dx=lamda/80; dt=dx; dlin=round(40.0/dx);  X_ros=-4; 

time=0; t=0; time_max=3000000; X_center=round(dlin/2.0)*dx;  h=1.0 ; m=1860.0; 


print("kolich shagov = ",dlin)
print("lamda = ",lamda)
print("dx,dt = ",dt)
print("k = ",kv)

plt.xlabel('Values')
plt.ylabel('Probability')


def f(x):
    x = i*dx - X_center
    if(x>=-dlin_p and x<=dlin_p):
        f=glub
    else:
            f=0.0
    return f 


def Psi(x):
    x = i*dx - X_center
    Psi=e**(-(x-X_ros)**2+kv*x*1j)                  
    return Psi 

U=[]

for i in range(dlin+1):                                                                                         
    x = i*dx - X_center                                                                                         
    U+=[(Psi(i*dx - X_center))]

r=[]

for i in range(1,dlin+1):                                                                                                                                                                         
    r+=[i*dx - X_center] 

l=0
for n in range(1,dlin-1):
    l+=(abs(U[n])**2 + abs(U[n+1])**2)/2



B=1.0+1j*dt/(2.0*m*dx**2) + 1j*f(x)*dt/2.0
C=A=1j*dt/(4.0*m*dx**2)

time1=199 
g=[]

pylab.ion()
#d=[]


while True: 

    time+=1 
    V=[]
    p=[]; p+=[0.0+0j];
    q=[]; q+=[0.0+0j];
    for i in range(1,dlin):
        x = i*dx - X_center
        D=(U[i] + (U[i+1]-2*U[i]+U[i-1])*1j*dt/(4*m*dx**2)-1j*dt*f(x)*U[i]/2)
        w2=B-C*p[i-1]
        p+=[A/w2]
        q+=[(D+A*q[i-1])/w2]
    V+=[0.0+0j]
    k=0.0+0j
    for i in range(dlin-1,0,-1):
        k=p[i]*k+q[i]
        V+=[k]
    V+=[0.0+0j]
    del U
    U=[]
    for i in range(dlin+1):
        U+=[V[dlin-i]] 
    del p; del q; del V;
    time1+=1
    if (time1==200):
        del g
        #del d
        g=[]
        l=0
        #d=[]
        for n in range(1,dlin-1):
            l+=(abs(U[n])**2 + abs(U[n+1])**2)/2
 
        for i in range(1,dlin+1):
            g+=[U[i].real]
            #d+=[abs(U[i])**2]
        time1=0
        pylab.clf()
        plt.plot([-dx*dlin/2,-dlin_p,-dlin_p,dlin_p,dlin_p,dx*dlin/2],[0,0,glub,glub,0,0])
        #kve=np.array(d)
        kx=np.array(r)
        ky=np.array(g)
        pylab.plot(kx,ky)
        #pylab.plot(kx,kve)
        plt.grid(True)
        plt.text(-9,1.2, 'k = ' );plt.text(-8,1.2, round(kv,2));plt.text(-5,1.2, 'psi = ' );plt.text(-3.5,1.2, round(l*dx,2));plt.text(0,1.2, '  lamda = ' );plt.text(3,1.2, lamda);plt.text(5,1.2, 't =' );plt.text(6.1,1.2, round(time*dt))
        plt.axis([-11,11, -1.1, 1.1])
        pylab.draw()
        fname = '_tmp%03d.png'%time
        plt.savefig(fname)

    if(time==time_max):                                        
        break
pylab.close()