# Applying neumann boundary conditions to diffusion equation solution in python [duplicate]

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
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()


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]