I have developed the code below to solve a non-linear scalar hyperbolic conservation law (Burgers' equation). I want to advance to solving a one dimensional system of Euler equations for Gas Dynamics (which is also a hyperbolic partial differential equation) in the similar manner. Can anyone advise on how to do this? Thanks.
""" """ from numpy import * from scitools.std import * import time import glob,os for filename in glob.glob('Burgers*.png'): os.remove(filename) N = 1001 M = 75; Dx = 1.5/M Dt = 0.0005 a = 1.0 def Initialise(): U = zeros(M) K = linspace(0.0,1.5, M+1) for i in range(M): if K[i] <= 0.5: U[i] = -0.5 elif K[i] >= 0.5 and K[i] <= 1 : U[i] = 1 elif K[i] >= 1: U[i] = 0 return U def Evolve_in_One_Timestep (U): #Create an array to store the new values of the solution U_new = zeros((M,), float32) for k in range (1, M-1): derivative = -0.25*(U[k+1]**2 - U[k-1]**2)*(Dt/Dx) + 0.5*sqrt(a)*(U[k+1] - 2*U[k] + U[k-1])*(Dt/Dx) U_new[k] = U[k] + derivative U_new = -0.5 U_new[M-1] = 0 return U_new #Initialising the data array U_data = Initialise() counter = 0 for n in range(0,N): U_data = Evolve_in_One_Timestep(U_data) t = n*Dt if (n%100==0): plotlabel= "t = " + str(n*Dt) #time labelling plot(linspace(0.0,0.75,M), U_data,'r0.3',label=plotlabel,legend='t=%4.2f' % t , hardcopy = 'Burgers%04d.png' % counter ) counter += 1 time.sleep(0.5) # a pause to control movie speed movie('Burgers*.png')