python编程练习：Engquist-Osher差分格式求解Burgers方程

一、题目

二、代码

from scipy.integrate import quad
import numpy as np


# 函数f(u)=1/2*u**2,故f`(u)=u
def f_positive(upp_value):
    # 积分f+(u)中需要使用的函数
    is_over_zero = int(upp_value > 0)
    return is_over_zero * upp_value


def f_negative(upp_value):
    # 积分f-(u)中需要使用的函数
    is_over_zero = int(upp_value > 0)
    return (1 - is_over_zero) * upp_value


def cal_fu(upp_value, form='+'):
    """
    :param upp_value:积分上限
    :param form: 可选值列表[‘+’，‘-’],决定返回f+(u)还是f-(u)
    :return: quad_r: f+(u)、f-(u)的值
    """
    quad_r = 0
    if form == "+":
        quad_r = quad(f_positive, 0, upp_value)[0]
    elif form == '-':
        quad_r = quad(f_negative, 0, upp_value)[0]
    return quad_r


def cal_next_step(pt, px):
    """
    该空间点在下一时间层的值
    :param pt: 当前点在网格的空间位置
    :param px: 当前点在网格的时间位置
    :param grid: 网格点上的值
    :return:
    """
    u_j_n = grid_value[pt][px]  # 计算U(j,n)
    u_jplus1_n = grid_value[pt][px+1]  # 计算U(j+1,n)
    u_jminus1_n = grid_value[pt][px-1]  # 计算U(j-1,n)
    temp_minus = cal_fu(u_jplus1_n, form='-') - cal_fu(u_j_n, form='-')
    temp_plus = cal_fu(u_j_n) + cal_fu(u_jminus1_n)
    u_j_nplus1 = u_j_n - grid_ratio * (temp_minus + temp_plus)
    grid_value[pt+1][px] = u_j_nplus1


if __name__ == "__main__":
    x_range = [-2, 2]  # 空间范围
    t_range = [0, 0.9]  # 时间范围
    delta_x = 0.1  # 空间步长
    delta_t = 0.01  # 时间步长
    grid_ratio = delta_t / delta_x  # 网格比
    grid_x = int((x_range[1] - x_range[0]) / delta_x) + 1  # 空间网格点数，此例中为41
    grid_t = int((t_range[1] - t_range[0]) / delta_t) + 1  # 时间网格点数，此例中为91

    # 考虑用列表grid_value来存储Ujn[[t=0.01],...,[t=0.9]]
    grid_value = np.zeros((grid_t, grid_x))  # 行代表某个时间、列代表某个空间

    # 将初始值t=0添加到grid_value中,即初始条件
    for i in range(grid_x):
        x_current = x_range[0] + delta_x * i
        if x_current > 0:
            grid_value[0][i] = 1
        else:
            grid_value[0][i] = -1

    # 将每一个时间层上的左右边界赋固定值
    grid_value[:, grid_x-1] = 1  # 右边界为1
    grid_value[:, 0] = -1  # 左边界为-1

    # 开始计算，时间上索引从0算到89，空间上索引从1算到39
    # 假设右边界必定收敛
    for i in range(0, grid_t-1):
        for j in range(1, grid_x-1):
            cal_next_step(i, j)
        grid_value[i, -1] = grid_value[i, -2]
    grid_value[-1, -1] = grid_value[-1, -2]

    # 仅将最后一个时间层的网格点数据保存到"2.txt"中
    np.savetxt('2.txt', grid_value[-1, :], fmt='%0.8f')

三、运行结果

绘图如下：