For the diffusion equation $$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial ^2 u(x,t)}{\partial x^2} + Cu(x,t) $$ with the boundary conditions $u(-\frac{L}{2},t)=u(\frac{L}{2},t)=0$ I've programmed the numerical solution into python correctly (I think).

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
    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.title('Numerical Solution of the Diffusion equation over time')

However now I have to change the right boundary condition into $u_x(\frac{L}{2})=0$ and I'm not really sure how to change my code to reflect this. If I do this the critical length should decrease and the function should start increasing exponentially, but everything I've tried usually does nothing to my plot - is there something wrong with my original code possibly? Any help is really appreciated, I've been trying for ages but can't seem to get it! Thanks!


If you approximate $u_x(\frac{L}{2},t) = 0$ with a second-order difference, you get:

$u_x(\frac{L}{2},t) \approx \frac{1}{2\Delta x}\left(u_s^k - u_{s-2}^k \right) = 0$,

which simplifies to $u_s^k = u_{s-2}^k$.

Then, write the equation for the last node, $u_{s-1}$, just like you would for an interior node:

$\frac{1}{\Delta t} \left(u_{s-1}^{k+1} - u_{s-1}^{k} \right) = \frac{D}{\Delta x^2} \left(u_s^k - 2 u_{s-1}^k + u_{s-2}^k \right) + C u_{s-1}^k$.

Note that the above expression involves $u_s$, which is outside the domain (a ghost node), so we need to eliminate it by using the Neumann BC equation obtained earlier.

Substituting $u_s^k = u_{s-2}^k$ into the equation for $u_{s-1}$ and simplifying gives:

$u_{s-1}^{k+1} = u_{s-1}^{k} + \frac{D \Delta t}{\Delta x^2} \left(-2 u_{s-1}^k + 2 u_{s-2}^k \right) + C \Delta t\ u_{s-1}^k$

So in your code, you need to replace the right Dirichlet BC with the following:

u[s-1,k+1]=(1+Beta-2*Alpha)*u[s-1,k] + 2*Alpha*u[s-2,k]
