I have been given a convection-diffusion ODE modeling the steady state temperature of a pipe (through which flows a fluid) as

$$-\frac{d}{dz}\left(\kappa \frac{dT}{dz} \right)+v\rho C\frac{dT}{dz}=Q(z)$$.

Here $T$ is the (steady state) temperature and it depends on the space variable $z$. All other parameters are constants and for the function $Q$ holds

$$Q(z)= \left\{\begin{matrix} 0,&\quad 0\leq z <a \\ Q_{0}\sin(\frac{(z-a)\pi}{b-a}),&\quad a\leq z \leq b \\ 0,&\quad b<z\leq L \end{matrix}\right.$$

As for the boundary conditions (BCs) we have a Dirichlet condition for the leftmost boundary and a Neumann condition for the rightmost boundary,

$$T(0)=T_{0}\quad \text{and}\quad -\kappa \frac{dT}{dz}(L)=\alpha(v)(T(L)-T_{out}),$$

where $T_{out}$ is the temperature outside the pipe to the right of $z=L$ and $\alpha$ is given by

$$\alpha(v) = \sqrt{\frac{v^{2}\rho^{2}C^{2}}{4}+\alpha_{0}^{2}}-\frac{v\rho C}{2}. $$

Now, I would like to discretize this ODE using a second order finite difference method (FDM), so I propose to use the central difference approximation which approximates the first and second derivatives as

$$\begin{align*} \frac{du}{dz}(z_{i})&=\frac{u(z_{i+1})-u(z_{i-1})}{2h} + O(h^{2})\\ \frac{d^{2}u}{dz^{2}}(z_{i})&=\frac{u(z_{i+1})-2u(z_{i})+u(z_{i-1})}{h^{2}} + O(h^{2}) \end{align*}$$

Finally, I would like to discretize the Neumann BC, but here things get tricky. If I use a central difference, I will obtain a point outside the domain (ghost point?). Could I use the condition that the temperature outside of the pipe at $z=L$ is $T_{out}$, or is this false? And also, how do I implement or force the Neumann condition to be fulfilled? Can I simply replace $dT/dz$ in the ODE with the BC when approximating the solution at $z=L$? This is what I do right now.

Interestingly enough, my solution seems to blow upp when I increase the number of steps, or mesh points, in the discretization of the $z$-domain. Look at these results:

enter image description here

enter image description here

The temperature completely sky rockets! Why is this? (If you would like me to post me code, I am more than willing to do so.)

Best regards,

  • 1
    $\begingroup$ I think that type right boundary condition is a Robin BC, instead of a Neumann one. How are you solving your system of equations? $\endgroup$
    – nicoguaro
    Commented Sep 15, 2019 at 23:37
  • $\begingroup$ Yes, that is correct, my fault. In Matlab, I allocate a zero matrix for all the equations I obtain from the central difference approximation. Since Q changes over time, so does the rhs of the equations. Essentially, I solve this system Ax=f, and with Matlabs backslash inverse function I get the solution as x=A\f (that is A^{-1}f). And as I stated already when I implement the Robin condition I simply replace the first derivative evaluated at z=L in the ODE, then I combine this with setting the function value of the point outside the domain to T_{out} (which doesnt feel right...) @nicoguaro $\endgroup$ Commented Sep 16, 2019 at 8:09

1 Answer 1


Just use the left-sided derivative at the right boundary. Then your Robin boundary condition at the last grid point $n$ is expressed as $ \kappa \left( \frac{1}{2} T_{n-2} −2 T_{n-1} + \frac{3}{2} T_n \right) = \alpha h (T_n - T_{out}) $

  • $\begingroup$ Great! So this will be the last equation of my system? Is this a second order approximation? $\endgroup$ Commented Sep 16, 2019 at 8:11

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.