Application to the tube chock problem

En el libro An Introducctión to Scientific Computing Twelve Computational Projects Solved with MATLAB  presentan doce proyectos, que si bien estan resueltos, son un buen ejercicio eso días en los que no hay deberes. Entre mis intereses están los métodos numéricos para resolver las ecuaciones que gobiernan la dinámica de los fluidos (CFD-Dinámica de fluidos Computacional). Los proyectos 10 y 12 tratan estas ecuaciones. Hay algunas cuestiones importantes en la solución numérica de estas ecuaciones, que tienen que ver con la matemática de las ecuaciones, soluciones discontinuas, condiciones iniciales discontinuas; ver la matemática necesaria para entender estas ecuaciones es una de mis expectativas del curso de EDP 2 que voy a llevar este ciclo. 

No tengo MATLAB por el momento, así que me las arregle, use SAGE y lo programe en Python; estoy obteniendo los resultados.

class Tubechock:
    def __init__(self):
        self.gamma=1.4
        self.cfl=0.95
        self.g1=line([])
        self.g2=line([])
        self.g3=line([])
        self.g4=line([])
        self.M=81
        self.x0=0.5
        self.rhoL=8.0
        self.rhoR=1.0
        self.pL=10.0/self.gamma
        self.pR=1.0/self.gamma
        self.UL=0
        self.UR=0
        self.rhoUL=self.rhoL*self.UL
        self.rhoUR=self.rhoR*self.UR
        self.EL=self.pL/(self.gamma-1)
        self.ER=self.pR/(self.gamma-1)
        self.dx=1/(self.M-1)
        self.dt=0
        self.x=[i*self.dx for i in range(self.M)]
        self.t=0
        self.w=[[0 for j in range(3)]for i in range(self.M)]
        self.F=[[0 for j in range(3)]for i in range(self.M)]
        self.FF=[[0 for j in range(3)]for i in range(self.M-1)]
        for i in range(self.M):
            if (self.x[i]<=self.x0):
                self.w[i]=[self.rhoL,self.rhoUL,self.EL]
            if (self.x[i]>self.x0):
                self.w[i]=[self.rhoR,self.rhoUR,self.ER]
    
    
    def plot_d(self):        
        xvsd=[[self.x[i],self.w[i][0]] for i in range(self.M)]
        self.g1=line(xvsd,rgbcolor=(3/4,1/2,5/8))
        self.g1.axes_labels(['$x$',r'$\rho$'])
        self.g1.fontsize(15)
        
    
    
    def plot_v(self):
        xvsv=[[self.x[i],self.w[i][1]/self.w[i][0]] for i in range(self.M)]        
        self.g2=line(xvsv,rgbcolor=(1/4,1/4,5/8))
        self.g2.axes_labels(['$x$','$v$'])
        self.g2.fontsize(15)
    
    
    def plot_p(self):
        xvsp=[]
        for i in range(self.M):
            k=(self.w[i][2]-0.5*self.w[i][1]*self.w[i][1]/self.w[i][0])*(self.gamma-1)
            xvsp.append([self.x[i],k])
        self.g3=line(xvsp,rgbcolor=(6/8,3/4,3/8))
        self.g3.axes_labels(['$x$','$p$'])
        self.g3.fontsize(15)
 
 
    def plot_T(self):
        xvsT=[]
        for i in range(self.M):
            k=(self.w[i][2]-0.5*self.w[i][1]*self.w[i][1]/self.w[i][0])*(self.gamma-1)
            xvsT.append([self.x[i],k/self.w[i][0]])
        self.g4=line(xvsT,rgbcolor=(3/4,3/5,1/3))
        self.g4.axes_labels(['$x$','$T$'])
        self.g4.fontsize(15)
        
    
    def graf(self):
        self.plot_d()
        self.plot_v()
        self.plot_p()
        self.plot_T()
        g=graphics_array([[self.g1,self.g2],[self.g3,self.g4]])
        g.show(frame=True, axes=True, figsize=[12,8])
        
    
    
    def Dt(self):
        uloc=[(self.w[i][1]/self.w[i][0]) for i in range(self.M)]
        press=[(self.gamma-1)*(self.w[i][2]-0.5*self.w[i][1]*uloc[i]) for i in range(self.M)]
        k=self.gamma*press[i]/self.w[i][0]        
        a=[((k)^(0.5)) for i in range(self.M)]
        max=[(abs(uloc[i])+a[i]) for i in range(self.M)]
        max.sort()
        self.dt=0.95*self.dx/max[self.M-1]
        
        
    
    def Fw(self):
        for i in range(self.M):            
            uloc  = self.w[i][1]/self.w[i][0]
            rhou2 = self.w[i][1]*uloc
            press = (0.4)*(self.w[i][2]-0.5*rhou2)
            self.F[i]=[self.w[i][1],rhou2+press,(self.w[i][2]+press)*uloc]
            
    def FFw(self):
        for i in range(self.M-1):
            pas=[0.5*(self.w[i][j]+self.w[i+1][j])-0.5*(self.dt/self.dx)*(self.F[i+1][j]-self.F[i][j]) for j in range(3)]
            uloc  = pas[1]/pas[0]
            rhou2 = pas[1]*uloc
            press = (0.4)*(pas[2]-0.5*rhou2)
            self.FF[i]=[pas[1],rhou2+press,(pas[2]+press)*uloc]
                
    
    
    def sol(self,tf=0.2):
        while (self.t<tf):                    
            self.Dt()
            self.Fw()
            self.FFw()
            for i in range(1,self.M-1):
                self.w[i][0]=self.w[i][0]-(self.dt/self.dx)*(self.FF[i][0]-self.FF[i-1][0])
                self.w[i][1]=self.w[i][1]-(self.dt/self.dx)*(self.FF[i][1]-self.FF[i-1][1])
                self.w[i][2]=self.w[i][2]-(self.dt/self.dx)*(self.FF[i][2]-self.FF[i-1][2])
            self.t=self.t+self.dt            
        
 
mysol=Tubechock()
mysol.sol()
mysol.graf()