计算科学 - 在python中将neumann边界条件应用于扩散方程解 - 吾爱随笔录

对于扩散方程

\frac{\partial u (x, t)}{\partial t} = D \frac{\partial^{2} u (x, t)}{\partial x^{2}} + C u (x, t)

$\frac{\partial u(x,t)}{\partial t} = D \frac{\partial ^2 u(x,t)}{\partial x^2} + Cu(x,t)$ 有边界条件

u (- \frac{L}{2}, t) = u (\frac{L}{2}, t) = 0

$u(-\frac{L}{2},t)=u(\frac{L}{2},t)=0$ 我已经正确地将数值解编程到 python 中（我认为）。

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)) #setting the function u
u[50,0]=1/ds # delta function
for k in range(0,t-1):
    u[0,k]=0 # boundary conditions
    u[s-1,k]=0
    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.title('Numerical Solution of the Diffusion equation over time')
plt.xlabel('x')
plt.ylabel('u(x,t)')
plt.show()

但是现在我必须将正确的边界条件更改为 $u_x(\frac{L}{2})=0$ 而且我不确定如何更改我的代码以反映这一点。如果我这样做，临界长度应该减少并且函数应该开始呈指数增长，但是我尝试过的一切通常对我的情节没有任何作用——我的原始代码可能有问题吗？任何帮助都非常感谢，我已经尝试了很长时间，但似乎无法得到它！谢谢！