I am considering the following diffusion equation:
$$\frac{\partial f}{\partial t} = \frac{\partial}{\partial x}[D(x,t)\frac{\partial f}{\partial x}]$$
over a grid x=[1,...,6]
with the following boundary conditions:
$$f|_{x=1} = 0 \quad \forall t$$
$$\partial f|_{x=6} = 0 \quad \forall t$$
To test for convergence, I am using an arbitrary solution (which maintains the boundary conditions):
$$f = \sin(\pi/10 * x - \pi/10)\sin(bt)$$
which after plugging into the diffusion equation introduces a source term which, when added to the original governing equation, yields a new equation for which the manufactured solution is an exact solution (easily found by substituting the artificial solution into the diffusion equation).
Running for $t=600$s, the analytic solution looks like the following:
But my numerical solution (using Crank-Nicolson) looks like the following, with an error at the Neumann boundary:
I believe my code is fine (listed below) and wondered if this could just be a case of my diffusion equation above not being well-posed for different boundary conditions on each end?
Both boundaries were formulated using ghost points, the right Neumann boundary specifically by:
$$\frac{f_{6+\Delta x} - f_{6-\Delta x}}{2\Delta x} = 0$$
Could this also be a cause of the error, and is there a better way to define a Neumann boundary?
Thank you for your help and I'm happy to upload some code if this would be helpful.
UPDATE: The error seems to occur due to incorrectly using the two-sides approximation for the Neumann right boundary. Adopting:
$$\frac{f_{6} - f_{6-\Delta x}}{\Delta x} = 0$$
instead, I now get what I want.
There is a still a small error at the boundary, but I believe this can be improved with an $\mathcal{O}(\Delta x ^2)$ approximation for the Neumann condition.
UPDATE 2: For reproducibility, the ungeneralised form of my diffusion equation is
$$\frac{\partial f}{\partial t} = x^2 \frac{\partial}{\partial x}[\frac{D(x,t)}{x^2}\frac{\partial f}{\partial x}]$$
and if we say $\bar{D} = \frac{D(x,t)}{x^2}$, I solve with Crank-Nicolson via the following:
$$\frac{f_{j}^{n+1} - f_{j}^{n}}{\Delta t} = \frac{x_j^2}{2}\bigg[\frac{\bar{D}_{j+\frac{1}{2}}^{n+\frac{1}{2}}(f_{j+1}^{n}-f_j^n) - \bar{D}_{j-\frac{1}{2}}^{n+\frac{1}{2}}(f_{j}^{n}-f_{j-1}^n)}{(\Delta x)^2} + \frac{\bar{D}_{j+\frac{1}{2}}^{n+\frac{1}{2}}(f_{j+1}^{n+1}-f_j^{n+1}) - \bar{D}_{j-\frac{1}{2}}^{n+\frac{1}{2}}(f_{j}^{n+1}-f_{j-1}^{n+1})}{(\Delta x)^2}\bigg]$$
My code for this scheme with the above boundary conditions is as follows:
def Crank_Nicolson(dt,nt,dL,L,f,Dlist,Q,):
'''
Must give D and Q for (nt+1) times
'''
T = f.copy()
lbc = T[0]
s = (0.5*dt/dL**2)
res = []
res.append(f)
for n in range(nt):
D = Dlist[n]
Dplus = Dlist[n+1]
Dl = np.array([L[i]**2 * 0.5*((D[i] + D[i-1])/2 + (Dplus[i] + Dplus[i-1])/2)/(L[i]-0.5*dL)**2 for i in range(1,len(L)-1)])
Dr = np.array([L[i]**2 * 0.5*((D[i] + D[i+1])/2 + (Dplus[i] + Dplus[i+1])/2)/(L[i]+0.5*dL)**2 for i in range(1,len(L)-1)])
Dc = np.array([x+y for x,y in zip(Dl,Dr)])
A = diags([-s*Dl[1:], 1+s*Dc, -s*Dr], [-1, 0, 1], shape=(len(L)-2, len(L)-2)).toarray()
B1 = diags([s*Dl[1:], 1-s*Dc, s*Dr], [-1, 0, 1], shape=(len(L)-2, len(L)-2)).toarray()
#Boundary conditions
b = np.zeros(len(L)-2)
b[0] = 2*s * 0.5*((D[1] + D[0])/2 + (Dplus[1] + Dplus[0])/2) \
/ (L[1]-0.5*dL)**2 \
*L[1]**2 * lbc
Dminus = 0.5*((D[-2] + D[-3])/2 + (Dplus[-2] + Dplus[-3])/2) \
/(L[-2]-0.5*dL)**2 *L[-2]**2
A[-1,-1] = 1 + (s* Dminus)
B1[-1,-1] = 1 - (s* Dminus)
A[-1,-2] = - s* Dminus
B1[-1,-2] = s* Dminus
Qnew = (Q[n]+Q[n+1])/2
Tn = T.copy()
# B = np.add(np.dot(Tn[1:-1],B1),b+ dt*(Qnew[1:-1])) is incorrect
B = np.add(np.dot(B1,Tn[1:-1]),b+ dt*(Qnew[1:-1]))
T[1:-1] = np.linalg.solve(A,B)
T[-1] = T[-2]
res.append(T.copy())
return T,dt, res, L, dL
UPDATE 2: Turns out the error at the right Neumann boundary was due to the code, in the line B = np.add(np.dot(Tn[1:-1],B1),b+ dt*(Qnew[1:-1]))
, since I forgot that MATRIX MULTIPLICATION IS NOT COMMUTABLE. Correcting this I finally get the desired result