Применение граничных условий Неймана к уравнению диффузии

Question

Я построил код для численного решения уравнения диффузии du / dt = D (d ^ 2 u / dx ^ 2) + Cu, где u - функция от x и t - я решил эточисленно и нанесено на график с помощью граничных условий диреклетов u (-L / 2, t) = u (L / 2, t) = 0, где критическая длина является значением до экспоненциального взрыва функции, к которому я разработалбыть пи.Я пытаюсь изменить правильное граничное условие на граничное условие Неймана 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()

Для дирихлета BC вы получаете: (вероятно, то же самое, что и раньше, поскольку для задачи диффузии это не будет иметь большого значения, но в любом случае было неверным),Затем вы меняете свое правое граничное условие для BC-Neumann BC

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

, поскольку я вижу, что вы используете центральную разностную схему, поэтому BC-Neumann BC утверждает, что du / dx = 0 на границе.Дискретизация этой производной с центральной разностной схемой на правой границе равна (u[s-1,k]-u[s-3,k])/dx = 0, поэтому u[s-1,k]=u[s-3,k].С этим изменением вы получите