离散Burgers方程使用阶跃函数时出现奇异矩阵错误求助

我在实现离散Burgers方程时遇到了一个棘手的问题：当使用阶跃函数作为初始条件时，程序会抛出singular matrix error，但换成平滑的初始条件（比如v = 0.4 + 0.2 * np.exp(-40 * (x - 0.5)**2)）时却能正常运行。阶跃函数理论上应该能得到合理结果，但我始终定位不到问题根源。

我的运行参数如下：

• 粘度值ν = 0.01
• 时间步长dt = 0.01
• 网格尺寸h = 0.01
• 最终时间T = 1

初始解u(x, 0)取x=0.1处的单位阶跃，边界条件设置为：u(0, t) = 1和u(1, t) = 0。

我是Python编程的新手，恳请各位帮忙排查问题。以下是我的实现代码：

离散Burgers方程实现 #

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

    # ua 为左边界条件
    # ub 为右边界条件 
    # nu 为运动粘度 
    m = ukp.size
    # f 用于存储空间中每个uk点的函数值
    f = np.zeros((m-2, 1))

    # 边界条件赋值
    uL = ua
    uR = ub

    # 左边界离散方程
    f[0] = (uk[0] - ukp[1])/dt + uk[0] * (uk[0] - uL)/h - nu * (uk[1] - 2*uk[0] + uL)/h**2

    # 内部节点的离散方程
    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

    # 右边界离散方程
    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

初始数据设置 #

def initialData(m):

    # 输入: 网格尺寸 m
    # 输出: 初始条件 u
    
    # 空间域端点
    xL = 0
    xR = 1
    
    h = (xR - xL) / (m-1)                      # 网格尺寸
    x = np.linspace(xL, xR, m).reshape((m, 1)) # 网格点
    v = np.zeros(len(x))

    v = my_step_function(x) # 初始数据
    
    # 边界条件提取
    uL = v[0]
    uR = v[m-1]
    
    return v

内容的提问来源于stack exchange，提问作者Nexi