Is the diffusion equation with Neumann and Dirichlet BCs well-posed?

Question

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

Answer 1

Yes, the problem with mixed boundary conditions is well posed. What's not clear to me is this:

• Why do you approximate the derivative via the two-sides approximation? Shouldn't it be enough to just the following? $$ \frac{f_6 - f_{6-\Delta x}}{\Delta x} = 0. $$

• In your animations, what is the size of $\Delta x$? Your curve looks very smooth, which can be either because you have a very small $\Delta x$, or because you take all of the data points $(x_i,f_i)$ and then put some curve through these points for the purposes of plotting. But that curve is not the data you computed and it might be that you are only seeing an artifact of plotting, rather than the data you're actually computing.

Answer 2

Problem well posed

Your problem is well posed.

On the discretization with Crank-Nicholson

I am not familiar with MMS, and I wonder how you got that ungeneralised form of the diffusion equation. Anyways, as I understand it, you are using the Crank-Nicholson method to discretize the following differential equation:

$\frac{\partial f}{\partial t} = x^2 \frac{\partial}{\partial x}[\frac{D(x,t)}{x^2}\frac{\partial f}{\partial x}]$

by means of the Crank-Nicholson method.

This equation is a coefficient variable diffusion equation, and I am not sure if C-N (Crank-Nicholson) is suited for that. I would first write the equation as a standard diffusion-advection equation, since Crank-Nicholson has been used many times for that purpose. As I'm not very sure of the above, I posted a question on this matter.

Using the product rule, the differential equation can be rewritten as a diffusion-advection equation:

$\frac{\partial f}{\partial t} = \frac{\partial}{\partial x}[D(x,t)\frac{\partial f}{\partial x}] - D(x,t)\frac{2}{x}\frac{\partial f}{\partial x}$

And this is the equation I would discretize with the C-N method. To that end, let us discretize the spatial domain (with central differences):

$\dot{f}_i = D_{i+\frac{1}{2}} \frac{f_{i+1}-f_i}{\Delta x^2} - D_{i-\frac{1}{2}} \frac{f_i-f_{i-1}}{\Delta x^2} - \frac{D_i}{x_i}\frac{f_{i+1}-f_{i-1}}{\Delta x} \equiv F_i$

where $\dot{f}$ refers to the derivative of $f$ with respect to time. Discretizing the time yields:

$\frac{f_i^{n+1}-f_i^n}{\Delta t} = \frac{1}{2} ( F_i^n + F_i^{n+1} )$

So, the fully discretized equations reads:

$\frac{f_i^{n+1}-f_i^n}{\Delta t} = \frac{1}{2} ( [D_{i+\frac{1}{2}}^n \frac{f_{i+1}^n-f_i^n}{\Delta x^2} - D_{i-\frac{1}{2}}^n \frac{f_i^n-f_{i-1}^n}{\Delta x^2} - \frac{D_i^n}{x_i}\frac{f_{i+1}^n-f_{i-1}^n}{\Delta x}] + [D_{i+\frac{1}{2}}^{n+1} \frac{f_{i+1}^{n+1}-f_i^{n+1}}{\Delta x^2} - D_{i-\frac{1}{2}}^{n+1} \frac{f_i^{n+1}-f_{i-1}^{n+1}}{\Delta x^2} - \frac{D_i^{n+1}}{x_i}\frac{f_{i+1}^{n+1}-f_{i-1}^{n+1}}{\Delta x}])$

On the Neumann boundary condition

The approximation proposed by @WolfgangBangerth is first order accurate (actually, in this particular case where the condition is set to 0, it is second order accurate, as pointed out by @whpowell96). However, your solution is second order accurate, and should be better.

How to impose the condition? We know that:

$\frac{f_{6+\Delta x} - f_{6-\Delta x}}{2\Delta x} = 0$

Therefore, you know that $f_{6+\Delta x} = f_{6-\Delta x}$. This is what you need to introduce in the last row of your matrix. In particular, recall equation:

$F_i = D_{i+\frac{1}{2}} \frac{f_{i+1}-f_i}{\Delta x^2} - D_{i-\frac{1}{2}} \frac{f_i-f_{i-1}}{\Delta x^2} - \frac{D_i}{x_i}\frac{f_{i+1}-f_{i-1}}{\Delta x}$

Let $j$ be such that $1 + j\Delta x =6$, and using the relation $f_{j + 1} = f_{j - 1}$ (from the Neuamnn boundary condition), you get:

$F_j = \frac{2}{\Delta x^2}D_{j-\frac{1}{2}} (f_{j-1}-f_j)$

where, in a similar way, we set $D_{j+\frac{1}{2}}$ equal to $D_{j-\frac{1}{2}}$.

UPDATE

According to this, both discretizations are $O(\Delta x^2)$, so both can be considered equivalent.