将纽曼边界条件应用于扩散方程

Question

我已经绘制了扩散方程 du/dt=D(d^2 u/dx^2) + Cu 的数值解的代码，其中 u 是 x 和 t 的函数 - 我已经解决了数值并用 direchtlet 边界条件 u(-L/2,t)=u(L/2,t)=0 绘制它，临界长度是函数以指数方式爆炸之前的值，我已经计算出成为圆周率。我正在尝试将正确的边界条件更改为 Neumann 边界条件 u_x(L/2)=0，这应该会减少临界长度并导致函数呈指数增长，但我不太确定该怎么做而且我的似乎不太正确-有人可以看看他们是否可以确定我哪里出错了吗？谢谢！

L=np.pi # value chosen for the critical length
s=101 # number of steps in x
t=10002 # number of timesteps
ds=L/(s-1) # step in x
dt=0.0001 # time step
D=1 # diffusion constant, set equal to 1
C=1 # creation rate of neutrons, set equal to 1
Alpha=(D*dt)/(ds*ds) # constant for diffusion term
Beta=C*dt # constant for u term

x = np.linspace(-L/2, 0, num=51)
x = np.concatenate([x, np.linspace(x[-1] - x[-2], L/2, num=50)]) # setting x in the specified interval

u=np.zeros(shape=(s,t))
u[50,0]=1/ds # delta function
for k in range(0,t-1):
    u[0,k]=0 #direchtlet boundary condition
    for i in range(1,s-1):
        u[i,k+1]=(1+Beta-2*Alpha)*u[i,k]+Alpha*u[i+1,k]+Alpha*u[i-1,k] # numerical solution 
    u[s-1,k+1]=u[s-2,k+1] # neumann boundary condition
    if k == 50 or k == 100 or k == 250 or k == 500 or k == 1000 or k == 10000: # plotting at times
        plt.plot(x,u[:,k])

plt.show()

Answer 1

首先，您需要在时间循环的开始和当前时间步中应用左右边界条件（而不是在k+1，就像您在右侧 BC 上所做的那样） .

import numpy as np
import matplotlib.pyplot as plt
L=np.pi # value chosen for the critical length
s=101 # number of steps in x
t=10002 # number of timesteps
ds=L/(s-1) # step in x
dt=0.0001 # time step
D=1 # diffusion constant, set equal to 1
C=1 # creation rate of neutrons, set equal to 1
Alpha=(D*dt)/(ds*ds) # constant for diffusion term
Beta=C*dt # constant for u term

x = np.linspace(-L/2, 0, num=51)
x = np.concatenate([x, np.linspace(x[-1] - x[-2], L/2, num=50)]) # setting x in the specified interval

u=np.zeros(shape=(s,t))
u[50,0]=1/ds # delta function
for k in range(0,t-1):
    u[0,k] = 0 # left direchtlet boundary condition
    u[s-1,k] = 0 # right dirichlet boundary condition

    for i in range(1,s-1):
        u[i,k+1]=(1+Beta-2*Alpha)*u[i,k]+Alpha*u[i+1,k]+Alpha*u[i-1,k] # numerical solution 
    if k == 50 or k == 100 or k == 250 or k == 500 or k == 1000 or k == 10000: # plotting at times
        plt.plot(x,u[:,k])

plt.savefig('test1.png')
plt.close()

对于 dirichlet BC，您会得到： （可能与以前相同，因为对于扩散问题，它不会有太大的不同，但无论如何都是不正确的）。然后你改变你的冯诺依曼BC的右边界条件

u[s-1,k] = u[s-3,k] # right von-neumann boundary condition

因为我看到您使用的是中心差分方案，所以 Von-Neumann BC 在边界处声明 du/dx=0。用右边界的中心差分方案离散这个导数是(u[s-1,k]-u[s-3,k])/dx = 0，所以u[s-1,k]=u[s-3,k]。有了这个改变，你得到