This question already has an answer here:

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!


marked as duplicate by Anton Menshov, nicoguaro Nov 28 '18 at 16:56

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.


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]

Not the answer you're looking for? Browse other questions tagged or ask your own question.