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:

enter image description here

But my numerical solution (using Crank-Nicolson) looks like the following, with an error at the Neumann boundary:

enter image description here

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.

enter image description here

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 = []
    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]


    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

enter image description here


2 Answers 2


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.

  • $\begingroup$ My $\Delta x =0.1$ and so is quite small. I was using the two-sides approximation following a web-page that I found, but now adopting your approach this is exactly where I am going wrong! I have added an update to the question adopting your solution $\endgroup$
    – R Thompson
    Commented Nov 25, 2019 at 15:51
  • 3
    $\begingroup$ There you go -- glad to hear! It might still be worth showing only the dots you compute, rather than some smooth interpolation of them. $\endgroup$ Commented Nov 25, 2019 at 21:58
  • $\begingroup$ It is worth noting that this approximation is only second order accurate because the Neumann boundary condition is $0$. If it were some constant, for instance, then it would only be first order accurate at the endpoint and a ghost point using central difference would be needed. Your numerical solution looks as if there is something off in either the coefficient matrix or something at the endpoint that is causing that effect $\endgroup$
    – whpowell96
    Commented Nov 26, 2019 at 20:29

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}}$.


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

  • $\begingroup$ Thank you for this detailed answer. I agree with your approach for turning into a diffusion-advection equation and then discretizing. The version I use is from a paper, which combined methods from both Press - Numerical recipes [e-maxx.ru/bookz/files/numerical_recipes.pdf] (page 1046) and [onlinelibrary.wiley.com/doi/abs/10.1002/cnm.879]. I was skeptical of it myself, but pressed on with the MMS. Since your method seems more robust, I will do a quick study to compare both forms of the C_N presented here and see what happens. This is actually very important for my research $\endgroup$
    – R Thompson
    Commented Nov 28, 2019 at 16:30
  • $\begingroup$ Did you impose the Neumann boundary as specified in the answer? $\endgroup$
    – mfnx
    Commented Nov 28, 2019 at 16:33
  • $\begingroup$ I will yes, I'm just doing some convergence testing with the first order accurate Neumann boundary and will then implement your second order accurate approach and do some more convergence testing $\endgroup$
    – R Thompson
    Commented Nov 28, 2019 at 16:35
  • $\begingroup$ @RThompson, if you obtain any result, let me know in the chat or something. It intrigues me. $\endgroup$
    – mfnx
    Commented Nov 28, 2019 at 20:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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