Численное решение нестационарного двумерного уравнения теплопроводности в python дает ошибку некорректно - PullRequest
Новая
       120

Численное решение нестационарного двумерного уравнения теплопроводности в python дает ошибку некорректно

0 голосов
/ 10 июля 2020

Я пытаюсь реализовать два численных решения. Прямой Эйлер и Рунге-Кутта второго порядка для нестационарного двумерного уравнения теплопроводности с периодическими c граничными условиями. Я использую трехточечную центральную разницу в обоих случаях для пространственной дискретизации.

Это мой код: 

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. Я довольно тщательно просмотрел код и не обнаружил явных опечаток или ошибок. Я чувствую, что это должно быть что-то незначительное, чего мне не хватает, но я не могу этого найти. Если вы случайно заметили ошибку или знаете что-то, я не могу помочь, или комментарий будет вам признателен.

Спасибо!

...