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

Question

我绘制了扩散方程的数值解的代码 du/dt=D(d^2 u/dx^2) + Cu 其中 u 是 x 和 t 的函数- 我已经用数值方法解决了它，并用 direchtlet 边界条件 u(-L/2,t)=u(L/2,t)=0 绘制了它，临界长度是函数之前的值呈指数级增长，我算出它是 pi。我正在尝试将正确的边界条件更改为 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]。通过此更改，您将获得

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

Applying neumann boundary conditions to the diffusion equation

python

numpy

matplotlib

math

numerical-methods