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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.