Upwind scheme

from scipy.integrate import solve_ivp

from numpy import zeros,exp,linspace,array,sin,cos,arange,pi

from numpy.linalg import norm

import matplotlib.pyplot as plt

# Decide scheme

scheme = "central"

# Define c

c = 1.5

# Define domain

delta_x     = 0.1

delta_t     = 0.05

final_t     = 3.0

x_values    = arange(0,4*pi+delta_x,delta_x)

t_values    = arange(0,final_t+delta_t,delta_t)

n = len(x_values)

m = len(t_values)

# Initialize empty vector to contain solution

phi_values = zeros((n,m))

# set initial condition

phi_values[:,0] = exp(-(x_values-pi)**2)

# advance in time

for j in range(1,m):

    # space sweep

    for i in range(n):

        if i == 0 or i == n-1:

            phi_values[i,j] = 0

        else:

            # for interior points

            phi_this        = phi_values[i,j-1]

            phi_east        = phi_values[i+1,j-1]

            phi_west        = phi_values[i-1,j-1]

            if scheme == "central":

                phi_values[i,j] = phi_this - \

                        (c*delta_t/(2*delta_x))*(phi_east-phi_west)

            elif scheme == "upwind":

                if c > 0:

                    phi_values[i,j] = 0

                elif c < 0:

                    phi_values[i,j] = 0

def plotstep(k):

    plt.plot(x_values[1:n-1],phi_values[1:n-1,k],label=f't = {k*delta_t:.2f}')

    return None

plt.figure(1)

plt.plot(x_values,phi_values[:,0],label=f't = 0')

plt.title(f'Scheme: {scheme}')

plotstep(30)

plotstep(50)

plt.legend()

plt.show()