Burger's equation (PDE) does not work with step function?

Question

I'm working on implementing the discretised Burger's equation. I am quite confused as to why it does not work with when using a step-function. When using a step-function, I am getting a singular matrix error but when using something else such as v = 0.4 + 0.2 * np.exp(-40 * (x - 0.5)**2), in the initial data function I do not.

The step function should provide results but I'm unable to identify where or why it's going wrong. I'm trying to run it with

• viscosity value ν = 0.01
• time-step t = 0.01
• grid size h = 0.01
• and final time T = 1.

The initial solution u(x, 0) is taken as a unit step at x = 0.1, boundary conditions as: u(0, t) = 1 and u(1, t) = 0.

Quite new to Python and notebook, would greatly appreciate insights into where I am going wrong.

I've provided a reproducible example: burgers.py

This is how I'm implementing the discretised Burger's equation:

def DiscretisedBurgers(uk, ukp, dt, h, nu, ua, ub):

    # ua is boundary condition
    # ub is boundary condition 
    # nu is kinematic viscosity 
    m = ukp.size
    # f to store values of the function for each point uk in space
    f = np.zeros((m-2, 1))

    # boundary conditions
    uL = ua
    uR = ub

    # left boundary
    f[0] = (uk[0] - ukp[1])/dt + uk[0] * (uk[0] - uL)/h - nu * (uk[1] - 2*uk[0] + uL)/h**2

    # Difference equations at each internal node
    for i in range(1, m-3):
        f[i] = (uk[i] - ukp[i+1])/dt + uk[i] * (uk[i] - uk[i-1])/h - nu * (uk[i+1] - 2*uk[i] + uk[i-1])/h**2

    # Difference equation at the right boundary
    f[m-3] = (uk[m-3] - ukp[m-2])/dt + uk[m-3] * (uk[m-3] - uk[m-4])/h - nu * (uR - 2*uk[m-3] + uk[m-4])/h**2

    return f

Initial data set up:

def initialData(m):

    # Input: grid size m
    # Output: initial condition u
    
    # spatial domain endpoints
    xL = 0
    xR = 1
    
    h = (xR - xL) / (m-1)                      # mesh size
    x = np.linspace(xL, xR, m).reshape((m, 1)) # grid points
    v = np.zeros(len(x))

    v = my_step_function(x) # initial data
    
    # Boundary conditions
    uL = v[0]
    uR = v[m-1]
    
    return v