I'm interested in solving the following 1D-advection equation using method of lines.

$$\frac{\partial C}{\partial t} = - v\frac{\partial C}{\partial x}$$ The spatial domain has been discretized into N nodes.Using forward difference for discretizing the spatial domain gives,

$$\frac{dC}{dt} = - v\frac{C_{i+1} - C_i}{\Delta x}$$

At the Nth node, Neumann boundary condition is applied

$$(\frac{dC}{dt})_N = 0$$

Implies, $C_{N+1} = C_N$

However, the following system doesn't converge. I am not sure if the above procedure is correct.

function convection()
global v;
N= 5;
C0 = [5 0 0 0 0];
v = 10000;
delx = 6;
tspan = [0 10];
[t C]  = ode15s(@(t,s) fun(t,s), tspan , C0);
function dC= fun(t,C)
    for i = 1:nnode-1
        dC(i,1) = -v*(C(i+1) - C(i))/delx; 
    dC(nnode,1) = 0;

I'd like to know whether boundary condition has to be applied on both ends. I have used Neumann BC only at the terminal node.


I think you are using a downwind- instead of an upwind finite difference. This leads to your code imposing a boundary condition where it is not allowed. The solution to your convection equation is basically (ignoring the left BC for the moment)

$$ C(x,t) = C_0(x - v t) $$

where $C_0$ is your initial value.

Thus, if $v > 0$, it is a rightward travelling signal. The value at the right end of your domain, in your case wherever your node $x_N$ is located, is completely determined from your initial value and left boundary condition. Imposing a Neumann BC there overdetermines the system and will lead to instability.

Second, your finite difference is biased in the wrong way. For positive $v$ (as in your example) it should be

$$ \frac{d C}{dt} = -v \frac{C_{i} - C_{i-1}}{\Delta x} $$

Changing this will also fix the issue with your boundary condition at $x_N$, since the upwind discretization won't need it.

  • $\begingroup$ Thanks a lot. I could impose Dirichlet BC at the left boundary by setting dC(1,1) = 0. I'd like to confirm whether Robin BC reduces to Dirichlet BC when there is no diffusion. $\endgroup$
    – Natasha
    Nov 19 '19 at 2:31
  • $\begingroup$ No, Robin BC is typically a combination of Dirichlet and Neumann BC. That is, something like u(0) = 2 u'(0), for example. I don't think diffusion has any effect here. $\endgroup$
    – Daniel
    Nov 19 '19 at 7:41
  • $\begingroup$ Sorry if I wasn't clear. When both convection and diffusion are involved, constant flux can be set at the left boundary using Robin BC, $J_L = vC_0 - DdC/dx$. However, when diffusive effects aren't involved I'm not sure how constant flux can be imposed at the left boundary. , Often constant flux BC is imposed using $J_L = vC_{at ghost node = -1}$ .This implies $vC_{-1} = vC_{0}$ in the absence of diffusion.Therefore, $C_{-1} =C_0$; dC/dt( at i=0) = 0. I am not sure if this is the right way to proceed if one has to implement a conatnt flux BC in the absence of diffusion. $\endgroup$
    – Natasha
    Nov 19 '19 at 10:25
  • 1
    $\begingroup$ Talking about fluxes suggests to me that you are now considering finite volumes, which are quite different to finite difference. There, you can impose a BC via a flux at the boundary. You will find more on this in lots of numerical maths textbooks. A personal favourite is the one by Dale Durran, but other tastes may vary. $\endgroup$
    – Daniel
    Nov 19 '19 at 11:53

Your Answer

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

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