report.py

# -*- coding: shift_jis -*-#------------------------------------------------------------------------------import numpy as npimport numpy.linalg as LAimport matplotlib.pyplot as pltimport os,sys#------------------------------------------------------------def title(name=None,make=0,suffix=''):    if name==None:        name = sys._getframe(1).f_code.co_name + suffix    #呼び出し元の関数名を取得する。    else:        name += suffix    print("="*10+name+"="*10)    if (make) and not os.path.isdir(name):        os.makedirs(name)    return name#------------------------------------------------------------getattrs=lambda obj,names:map(lambda name:getattr(obj,name),names.split())setattrs=lambda obj,names,attrs:[setattr(obj,name,attrs[i]) for i,name in enumerate(names.split())]#------------------------------------------------------------------------------class FDM:    '''移流方程式のEuler, RK2,RK4による時間発展シミュレーション空間方向は常に中央差分を使用'''#------------------------------------------------------------------------------    def __init__(self,**kargs):        '''FDM Objectの初期設定'''        scheme=kargs.get('scheme','Euler')  #EulerとFTCSは同じである。        assert scheme in ['Euler','RK2','RK4'], scheme        A=kargs.get('A',1.0) #無次元化された移流係数(Courant数)  (dt=dx=1と考える)        N = kargs.get('N',256)        psi=np.zeros((N,),float) #振幅用buffer        f=lambda psi1: 0.5*A*(np.roll(psi1,1)-np.roll(psi1,-1))  #(S-1/S)/2  (空間は常に中央差分)        setattrs(self,'scheme N psi f',[scheme,N,psi,f])#------------------------------------------------------------------------------    def set_psi(self):        '''初期振幅を設定する。'''        X=np.linspace(0,1,self.N)        self.psi[:]=np.exp(-((X-0.5)*10.0)**2)#------------------------------------------------------------------------------    def evolve(self):        '''1-step分時間発展させる。'''        scheme,N,psi,f=getattrs(self,'scheme N psi f')        if scheme=='Euler':            psi+=f(psi)        elif scheme=='RK2':            dpsi1=f(psi)            dpsi2=f(psi+dpsi1)            psi+=(dpsi1+dpsi2)/2        elif scheme=='RK4':            ############### ここにRK4によるpsiの更新式を記述 ###############            dpsi            pass#------------------------------------------------------------------------------def test_logistic(**kargs):    '''Logistic方程式の差分解法'''    suffix=kargs.get('suffix','')    folder=title(suffix=suffix,make=1)    x=kargs.get('x',0.001)  #初期値    dt=kargs.get('dt',1)    tmax=kargs.get('tmax',20)    Nplot=kargs.get('Nplot',2)    f=lambda x: x*(1-x)    xEuler=x   #xEuler: Euler用、x: RK4用    T,XEuler,X=[],[],[]  #T: 時刻用, XEuler: Euler用, X: RK4用    F1=plt.figure()    A1=F1.add_subplot(111)    for t in range(tmax+1):        T.append(t*dt)        XEuler.append(xEuler)        X.append(x)        if t%Nplot==0:            print(t)            Xtheory= np.exp(T) / (np.exp(T) + 999)            #######ここで、Xtheory(理論解)を設定 #############            A1.plot(T,Xtheory,T,XEuler,T,X)            A1.grid()            A1.set_xlim(T[0],T[-1]+1)            A1.set_ylim(-0.1,1.1)            A1.legend('theory Euler RK4'.split())            F1.savefig(os.path.join(folder,'x_%07d.png' % t))            A1.cla()        xEuler+=f(xEuler)*dt        ############### ここにRK4によるxの更新式を記述 ###############        dx1=f(x)*dt        dx2=f(x+dx1/2)*dt        dx3=f(x+dx2/2)*dt        dx4=f(x+dx3)*dt        x += (dx1+dx2*2+dx3*2+dx4)/6#------------------------------------------------------------------------------def test_FDM_multi(**kargs):    suffix=kargs.get('suffix','')    folder=title(suffix=suffix,make=1)    N= kargs.setdefault('N',256)    tmax=kargs.get('tmax',10000)    Nplot=kargs.get('Nplot',20)    y_range = kargs.get('y_range',[-0.5,1.5])    tit=kargs.get('tit','')    legend=kargs.get('legend',[])#...........................    FDMs=[]    multi=kargs.get('multi')    if multi:        keys, values_list= multi        if not isinstance(keys,list):            keys, values_list= [keys],[[values] for values in values_list]            #listでない場合はlist化する。        for values in values_list:            kargs.update(dict(zip(keys, values)))            FDMs.append(FDM(**kargs))    else:        FDMs.append(FDM(**kargs))    M=len(FDMs)    for i in range(M):        FDMs[i].set_psi()#...........................    F1=plt.figure()    A1=F1.add_subplot(111)    X=range(N)    for t in range(tmax+1):        if t % Nplot==0:            print(t)            Y=np.vstack([FDMs[i].psi for i in range(M)])            A1.plot(X,Y.T)            A1.grid()            A1.legend(legend)            A1.set_xlim(0,N)            A1.set_ylim(*y_range)            A1.set_title(tit+'  [t=%d]' % t)            F1.savefig(os.path.join(folder,'psi_%07d.png' % t))            A1.cla()        for i in range(M):            FDMs[i].evolve()#------------------------------------------------------------------------------def test_stability_multi(**kargs):    suffix=kargs.get('suffix','')    folder=title(suffix=suffix,make=1)    N=kargs.get('N',180) #plotする点の数    t_range=np.arange(0.0,np.pi*2,np.pi/N)    name=kargs.get('name','nomame')  #PNGファイル名    legend=kargs.get('legend',[])    tit=kargs.get('tit','')#.........................    kargs_list=[None]  #最初に単位円を描くのに対応して、１個dummyを入れておく。    multi=kargs.get('multi')    if multi:        keys, values_list= multi        if not isinstance(keys,list):            keys, values_list= [keys],[[values] for values in values_list]            #listでない場合はlist化する。        for values in values_list:            kargs.update(dict(zip(keys, values)))            kargs_list.append(dict(kargs))    else:        kargs_list.append(dict(kargs))#...............................    colors=['y']+['r','g','b','m','c']*10    F1=plt.figure()    A1=F1.add_subplot(111)    for i,kargs in enumerate(kargs_list):        if kargs is None:            G = lambda s:s            scheme = "円"            A = None        else:            scheme=kargs.get('scheme','Euler')            assert scheme in ['Euler','RK2','RK4'], scheme            A=kargs.get('A',0.5) #無次元化された移流係数            F=lambda s: -A*(s-1/s)/2#..................................            if scheme=='Euler':                GF = lambda F: 1+F                G=lambda s: GF(F(s))            elif scheme=='RK2':                GF= lambda F: 1+F+F**2/2                G = lambda s: GF(F(s))            elif scheme=='RK4':                GF= lambda F: 1+F+F**2/2+F**3/6+F**4/24                G = lambda s: GF(F(s))        X,Y=[],[]        check = False        for t in t_range:            s1 = np.exp(1.0j*t)            s2 = G(s1)            X.append(s2.real)            if abs(s2.real**2 + s2.imag**2) > 1.001:                # print(s2, s2.real**2 + s2.imag**2)                check = True        if not check:            print(scheme, A, "OK!")        else:            print(scheme, A, "無理")            # X.append(s2.real)            # Y.append(s2.imag)        # A1.scatter(X,Y,s=1,color=colors[i])#..............................    # A1.grid()    # A1.set_title('stability check:'+tit)    # A1.legend(['unit circle']+list(legend))    #legendがmap objectで与えられる場合も許容するためlistにcast    # F1.savefig(os.path.join(folder,name+'.png'))    # A1.cla()   #--------------------------------------------------------------------------------test_logistic(x=0.001)multi=['scheme','Euler RK2 RK4'.split()]#test_FDM_multi(multi=multi,legend=multi[1],tmax=10000,A=1.0,suffix='_A1_scheme@')multi,legend=['A',[1.0, 2.0]], '1.0 2.0'.split()#test_FDM_multi(multi=multi, legend=legend,scheme='RK4',tmax=10000,Nplot=100,dt=1.0,suffix='_RK4_A@')multi=['scheme',['Euler','RK2','RK4']]for A in range(270, 300):    A *= 0.01    print(A, " : ")    test_stability_multi(multi=multi,legend=multi[1],A=A, tit='A=%f' % A, name='scheme@_A%d' % int(A*100))