溢流：（34，＆＃x27;结果太大了）当数值求解非线性PDE时

Question

修复：问题是CFL条件不够低。

我应该通过将其应用于SOD的冲击管问题来比较不同的数值方法。当我使用输入参数el = 1且er = 0.1时（对应于左侧的0.4的压力为0.4，右侧为0.04），它起作用，但是其他所有人都使用1和0.1的压力，这意味着EL = 2.5，ER = 0.25 。但是，当我使用这些值时，我会遇到溢出错误，该如何解决？这是代码：

import numpy as np
import matplotlib.pyplot as plt


def init_func(domrange, ul, ur, rhol, rhor, El, Er):
    x = [[rhol, rhol*ul, El] if i < domrange/2 else [rhor, rhor*ur, Er] for i in range(domrange)]
    return x

def func(rho, mom, ener):
    gamma = 1.4
    p = (gamma-1)*(ener - (1/2)*(mom**2)/rho)
    f1 = mom
    f2 = (mom**2)/rho + p
    f3 = (mom/rho)*(ener + p)
    return [f1, f2, f3]

def g(l1, l2, i, deltax, deltat): # test different numerical methods, l1 = [density, momentum, energy], g(l1[i], l2[i]) = 1/2(f(l1)[i] + f(l2)[i]) - deltax/(2*deltat) (l2[i] - l1[i])
    val = (1/2)*(func(l1[0], l1[1], l1[2])[i] + func(l2[0], l2[1], l2[2])[i]) - (deltax/(deltat*2))*(l2[i] - l1[i])
    return val

def numsim(cfl, domain, deltax, alpha, ul, ur, rhol, rhor, El, Er): # creates matrix of data points along the x axis for a number of timesteps
    time = 2
    deltat = cfl * deltax / alpha
    timesteps = int(time / deltat)
    domrange = int(domain/deltax)
    func_t = [init_func(domrange, ul, ur, rhol, rhor, El, Er)] # [[[density, momentum, energy], ...]]

    for t in range(1, timesteps):
        func_x = []
        for x in range(domrange):
            func_i = []
            for i in range(3):
                func_p1 = func_t[t-1][x % domrange][i] - alpha*deltat/deltax * (g(func_t[t-1][x % domrange], func_t[t-1][(x+1) % domrange], i, deltax, deltat) - g(func_t[t-1][(x-1) % domrange], func_t[t-1][x % domrange], i, deltax, deltat))
                func_i.append(func_p1)
            func_x.append(func_i)
        func_t.append(func_x)

    return func_t

def uideo(matrix, deltax): # plot the graphs in sequence with a timer
    domrange = len(matrix[0])
    start = int((domrange*deltax/2 - 5)/deltax)
    stop = int((domrange*deltax/2 + 5)/deltax)
    xas = np.subtract(np.multiply(range(start, stop), deltax), domrange*deltax/2)
    N = len(matrix)
    
    for i in range(0, N, int(1/deltax/10)):
        rho = [matrix[i][x][0] for x in range(start, stop)]
        plt.clf()
        plt.cla()
        plt.plot(xas, rho)
        plt.pause(0.01)
    plt.show()

def plots(matrix, deltax):
    domrange = len(matrix[0])
    gamma = 1.4
    start = int((domrange*deltax/2 - 5)/deltax)
    stop = int((domrange*deltax/2 + 5)/deltax)
    xas = np.subtract(np.multiply(range(start, stop), deltax), domrange*deltax/2)
    N = len(matrix)

    rho = [matrix[N-1][x][0] for x in range(start, stop)]
    u = np.divide([matrix[N-1][x][1] for x in range(start, stop)], [matrix[N-1][x][0] for x in range(start, stop)]) # u = (mom/rho)
    E = [matrix[N-1][x][2] for x in range(start, stop)]
    p = (gamma - 1)*(np.subtract(E, (1/2)*np.divide((np.square(np.multiply(rho, u))), rho))) # p = (gamma - 1)(E - 1/2(mom)**2/rho)

    fig, axs = plt.subplots(2, 2)
    axs[0, 0].plot(xas, rho)
    axs[0, 0].set_title('rho')
    axs[0, 1].plot(xas, u, 'tab:orange')
    axs[0, 1].set_title('u')
    axs[1, 0].plot(xas, E, 'tab:green')
    axs[1, 0].set_title('E')
    axs[1, 1].plot(xas, p, 'tab:red')
    axs[1, 1].set_title('p')
    fig.tight_layout()

    for ax in axs.flat:
        ax.set(xlabel='x')

    # Hide x labels and tick labels for top plots and y ticks for right plots.
    #for ax in axs.flat:
    #    ax.label_outer()
    plt.show()

#uideo(numsim(cfl=0.5, domain=40.0, deltax=0.01, alpha=1.0, ul=0.0, ur=0.0, rhol=1.0, rhor=0.125, El=1.0, Er=0.1), deltax=0.01)
plots(numsim(cfl=0.5, domain=40.0, deltax=0.01, alpha=1.0, ul=0.0, ur=0.0, rhol=1.0, rhor=0.125, El=2.5, Er=0.25), deltax=0.01)

这是错误： Trackback（最近的最新电话）： 第11行，在Func中 p =（gamma -1）（enery - （1/2）（MOM ** 2）/rho） 溢流：（34，“结果太大”）