python - python中非定常二维热方程的数值解不正确地产生错误

Question

我正在尝试实现两个数值解决方案。具有周期性边界条件的非定常二维热方程的正向 Euler 和二阶 Runge-Kutta。在这两种情况下，我都使用 3 点中心差进行空间离散化。

这是我的代码：

def analytical(npx,npy,t,alpha=0.1): #Function to create analytical solution data
    X = numpy.linspace(0,1,npx,endpoint=False) #X spatial range
    Y = numpy.linspace(0,1,npy,endpoint=False) #Y Spatial Range
    uShape = (1,numpy.shape(X)[0],numpy.shape(Y)[0]) #Shape of data array
    u = numpy.zeros(uShape) #Allocate data array
    m = 2 #m and n = 2 makes the function periodic in 0->1 
    n = 2
    for i,x in enumerate(X): #Looping through x and y to produce analytical solution
        for j,y in enumerate(Y):
            u[0,i,j] = numpy.sin(m*pi*x)*numpy.sin(n*pi*y)*numpy.exp(-(m*m+n*n)*pi*pi*alpha*t)
    return u,X,Y

def numericalHeatComparisonFE(): #Numerical function for forward euler
    arraysize = 10 #Size of simulation array in x and y
    d = 0.1 #value of pde coefficient alpha*dt/dx/dx
    alpha = 1 #thermal diffusivity
    dx = 1/arraysize #getting spatial step
    dt = float(d*dx**2/alpha) #getting time step
    T,x,y = analytical(arraysize,arraysize,0,alpha=alpha) #get analytical solution
    numerical = numpy.zeros((2,)+T.shape) #create numerical data array
    ns = numerical.shape #shape of numerical array aliased
    numerical[0,:,:,:] = T[:,:,:] # assign initial conditions to first element in numerical
    error = [] #create empty error list for absolute error
    courant = alpha*dt/dx/dx #technically not the courant number but the coefficient for PDE
    for i in range(1,20):#looping through twenty times for testing - solving FE each step
        T,x,y = analytical(arraysize,arraysize,i*dt,alpha=alpha)
        for idx,idy in numpy.ndindex(ns[2:]):
            dxx = numerical[0,0,idx-1,idy]+numerical[0,0,(idx+1)%ns[-2],idy]-2*numerical[0,0,idx,idy] #X direction diffusion
            dyy = numerical[0,0,idx,idy-1]+numerical[0,0,idx,(idy+1)%ns[-1]]-2*numerical[0,0,idx,idy] #Y direction diffusion
            numerical[1,0,idx,idy] = courant*(dxx+dyy)+numerical[0,0,idx,idy] #Update formula
        error.append(numpy.amax(numpy.absolute(numerical[1,:,:,:]-T[:,:,:])))#Add max error to error list
        numerical[0,:,:,:] = numerical[1,:,:,:] #Update initial condition
    print(numpy.amax(error))

def numericalHeatComparisonRK2():
    arraysize = 10 #Size of simulation array in x and y
    d = 0.1 #value of pde coefficient alpha*dt/dx/dx
    alpha = 1 #thermal diffusivity
    dx = 1/arraysize #getting spatial step
    dt = float(d*dx**2/alpha) #getting time step
    T,x,y = analytical(arraysize,arraysize,0,alpha=alpha) #get analytical solution
    numerical = numpy.zeros((3,)+T.shape) #create numerical data array
    ns = numerical.shape #shape of numerical array aliased
    numerical[0,:,:,:] = T[:,:,:] # assign initial conditions to first element in numerical
    error = [] #create empty error list for absolute error
    courant = alpha*dt/dx/dx #technically not the courant number but the coefficient for PDE
    for i in range(1,20): #Test twenty time steps -RK2
        T,x,y = analytical(arraysize,arraysize,i*dt,alpha=alpha)
        for idx,idy in numpy.ndindex(ns[2:]): #Intermediate step looping through indices
            #Intermediate
            dxx = numerical[0,0,idx-1,idy]+numerical[0,0,(idx+1)%ns[-2],idy]-2*numerical[0,0,idx,idy]
            dyy = numerical[0,0,idx,idy-1]+numerical[0,0,idx,(idy+1)%ns[-1]]-2*numerical[0,0,idx,idy]
            numerical[1,0,idx,idy] = 0.5*courant*(dxx+dyy)+numerical[0,0,idx,idy]
            
        for idx,idy in numpy.ndindex(ns[2:]): #Update step looping through indices
            #RK Step
            dxx = numerical[1,0,idx-1,idy]+numerical[1,0,(idx+1)%ns[-2],idy]-2*numerical[1,0,idx,idy]
            dyy = numerical[1,0,idx,idy-1]+numerical[1,0,idx,(idy+1)%ns[-1]]-2*numerical[1,0,idx,idy]
            numerical[2,0,idx,idy] = courant*(dxx+dyy)+numerical[0,0,idx,idy]
        error.append(numpy.amax(numpy.absolute(numerical[2,:,:,:]-T[:,:,:]))) #Add maximum error to list
        numerical[0,:,:,:] = numerical[2,:,:,:] #Update initial conditions
    print(numpy.amax(error))

if __name__ == "__main__":
    numericalHeatComparisonFE()
    numericalHeatComparisonRK2()

运行代码时，我希望 RK2 的最大错误应该小于 FE，但我得到 0.0021498590913591187 了 FE 和 0.011325197051528346 RK2。我已经非常彻底地搜索了代码，没有发现任何明显的拼写错误或错误。我觉得它必须是我想念的小东西，但我似乎找不到它。如果您碰巧发现错误或知道我没有帮助的内容或发表评论，我们将不胜感激。

谢谢！

标签： pythonmathphysicsnumerical-methodsdifferential-equations