In [ ]:

    

"""
Copyright (c) 2014 Kunal Tyagi
GNU GPLv3
"""


    

In [19]:

    

"""
Program to solve
u_t + a*u_x = 0
u(t+1,x)-u(t,x) = -(a*delta_t/delta_x)*(u(t,x) - u(t,x-1))
"""
from numpy import *
import pylab


    

In [20]:

    

# initialisation of constants
LEN_X = 1.0
LEN_T = 1.0

DELTA_X = 0.01
DELTA_T = 0.01
CFL = 0.5

# derived constants
a = CFL*DELTA_X/DELTA_T
num_x = int(LEN_X/DELTA_X)
num_t = int(LEN_T/DELTA_T)


    

In [21]:

    

# matrix
u = zeros([num_t, num_x])
f = zeros([num_t, num_x])


    

In [52]:

    

# t=0 initialisation
for i in range(0, num_x):
    u[0][i] = 0
    f[0][i] = 0
    if i > 20 and i < 30:
        u[0][i] = 1 #sin(pi*(i-20)/10)


    

In [53]:

    

for i in range(0, num_t-1):
    for j in range(0, num_x):
        # u[i+1][j] = (1+CFL)*u[i][j] + CFL*u[i][j-1]
        # forward difference is unstable
        u[i+1][j] = u[i][j] - CFL*(u[i][j] - u[i][j-1])
plot(u[99], 'b')
plot(u[49], 'r')
plot(u[0], 'g')


    

    
Out[53]:





[<matplotlib.lines.Line2D at 0x7f4ec8978cd0>]






    

In [54]:

    

for i in range(0, num_t-1):
    for j in range(1, num_x-1):
        # Lax Wendroff
        u[i+1][j] = u[i][j] - CFL/2*(u[i][j+1] - u[i][j-1]) + (CFL*CFL)/2*(u[i][j+1] - 2*u[i][j] + u[i][j-1])
plot(u[99], 'b')
plot(u[49], 'r')
plot(u[0], 'g')


    

    
Out[54]:





[<matplotlib.lines.Line2D at 0x7f4ec8929fd0>]






    

In [55]:

    

for i in range(0, num_t-1):
    for j in range(2, num_x-1):
        # Beam Warming
        u[i+1][j] = u[i][j] - CFL*(u[i][j] - u[i][j-1]) - (CFL*(1-CFL))/2*(u[i][j] - 2*u[i][j-1] + u[i][j-2])
plot(u[99], 'b')
plot(u[49], 'r')
plot(u[0], 'g')


    

    
Out[55]:





[<matplotlib.lines.Line2D at 0x7f4ec87ebc10>]






    

In [56]:

    

b = zeros([num_t, num_x])
l = zeros([num_t, num_x])
for i in range(0, num_t-1):
    for j in range(2, num_x-1):
        # FROMM
        l[i+1][j] = u[i][j] - CFL/2*(u[i][j+1] - u[i][j-1]) + (CFL*CFL)/2*(u[i][j+1] - 2*u[i][j] + u[i][j-1])
        b[i+1][j] = u[i][j] - CFL*(u[i][j] - u[i][j-1]) - (CFL*(1-CFL))/2*(u[i][j] - 2*u[i][j-1] + u[i][j-2])
        u[i+1][j] = (b[i+1][j] + l[i+1][j])/2
plot(u[99], 'b')
plot(u[49], 'r')
plot(u[0], 'g')


    

    
Out[56]:





[<matplotlib.lines.Line2D at 0x7f4ec8743c50>]






    

In [ ]: