forward backward central differencing approximation

Derivatives applied to thermodynamic properties

When dealing with engineering problems we often face (partial) derivatives. Since most functions can not be differenciated analytically, we need numerical methods. This blog post deals with the numerical differentiation and its application to thermodynamic properties.

Theoretical Basics and Numerical error

First, let us start with a bit of theory. If you are only interested in the application, just continue reading. In case you wanna have a look intothe theory, hit the Read More button.

First of all, the derivative of a function f(x) can be defined as:

(1)   \begin{equation*} f'(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h}, \end{equation*}

Let us have a look at the Taylor series to understand this context. A Taylor series represents a function by an infinite sum of terms. So, lets expand the Taylor row for a function f at x:

(2)   \begin{equation*} f(x + h) = f(x) + h \cdot f'(x) + \frac{h^2}{2} \cdot f''(x) + \frac{h^3}{6} \cdot f'''(\xi) + \cdots. \end{equation*}

Here, \xi is in the interval (x, x+h). Solving for f'(x) leads to

(3)   \begin{equation*} f'(x) = \frac{f(x+h)-f(x)}{h}-\frac{h}{2} \cdot f''(x) -\mathcal{O}(h)^2. \end{equation*}

Here \mathcal{O} stands for the remainder of the higher order terms. Hence, the first error term is proportional to h. This means, the forward difference method is of first order. Additionally, another first order method is the backward difference method. It is not explained here in detail.

Let us consider a second Taylor series:

(4)   \begin{equation*} f(x - h) = f(x) - h \cdot f'(x) + \frac{h^2}{2} \cdot f''(x) - \frac{h^3}{6} \cdot f'''(\xi) + \cdots. \end{equation*}

By subtracting the first Taylor series from the second we get:

(5)   \begin{equation*} f(x + h) - f(a - h)= 2h \cdot f'(x) + \frac{h^3}{3} \cdot f'''(x) + \mathcal{O}(h)^5. \end{equation*}

Afterwards, let us solve that for f'(x)

(6)   \begin{equation*} f'(x) = \frac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{6} \cdot f'''(x) + \mathcal{O}(h)^4. \end{equation*}

Since the first error term is proportional to h^2, the central difference method is of second order.

In an exact differentiation, h would approach 0. While differentiating numerically, h is finite. The truncation error is the difference between a derivative and its finite difference representation.

Furthermore, we need to consider the roundoff error. The roundoff error represents the difference between an exact number and its numerical representation. For example, a variable with the data type double can hold approximately 16 significant figures. Simply put, it arises by chopping off numbers. The smaller h the larger the impact of the roundoff error.

Finite difference representation of derivatives

In order to calculate the derivative of a function numerically, we can use its finite difference representation. The forward difference approximation is of first order and can be defined as:

(7)   \begin{equation*} f'(x) \approx \frac{f(x+h)-f(x)}{h}. \end{equation*}

Furthermore, the backward difference approximation is also of first order and can be defined as:

(8)   \begin{equation*} f'(x) \approx \frac{f(x)-f(x-h)}{h}. \end{equation*}

Additionally, the central difference approximation is of second order and can be defined as:

(9)   \begin{equation*} f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}. \end{equation*}

Differencing approximations

Specific heat capacity

The specific heat capacity at constant pressure is defined as

(10)   \begin{equation*} c_p = \left. \frac{\partial h }{\partial T} \right|_p \end{equation*}

This means, the enthalpy h is derivated with respect to the temperature T under constant pressure p.

CoolProp is a fantastic open source library to receive thermodynamic properties. I use it extensively and there is a lot of functionality. Since it is even possible to calculate derivatives, we use it for benchmarking reasons.

Applying central difference to specific heat

We can calculate the specific heat capacity c_p by applying the central difference approximation as follows:

(11)   \begin{equation*} c_p = \frac{(h+\Delta h)-(h-\Delta h)}{T(p, h+\Delta h) - T(p, h-\Delta h)}. \end{equation*}

Hence, we need to be able to calculate a temperature T from the enthalpy h and pressure p (T = f(p, h)). We use a CoolProp function in order to save us some time.

Furthermore, it is super important to choose an adequate \Delta h by which the enthalpy is pertubed. Therefore, we use the machine epsilon \epsilon_m and define \Delta h as

(12)   \begin{equation*} \Delta h = \sqrt[p+1]{\epsilon_m} \cdot h. \end{equation*}

Since we use the central difference approximation our order is 2 (p = 2).

Writing the Python script

So let us go over the Python script.

  • Import the dependencies CoolProp and numpy
  • Input parameter: Nitrogen, 100 bar, 300 K
  • Get the machine epsilon
  • Specify the function specificHeat
  • Call the own and the CoolProp function
#Author: johannes <>, 14.12.17
#License: MIT License (

import numpy as np
import CoolProp.CoolProp as CP

fluid = "Nitrogen"

#Machine epsilon
eps = np.finfo(float).eps

def specificHeat(h, p):
    pert = pow(eps, 0.33333)*h
    num = (h + pert)-(h - pert)
    denum = CP.PropsSI("T", "P", p, "H", h+pert, fluid)-CP.PropsSI("T", "P", p, "H", h-pert, fluid)
    return num/denum

#Input data
p = 100e5 #Pa
T = 300  #K
h = CP.PropsSI("H", "P", p, "T", T, fluid) #J/kg
#Function calls  
cp = specificHeat(h, p)      
cp_DerivativeCP = CP.PropsSI('d(H)/d(T)|P','P',p,'T',T,fluid)

print cp
print cp_DerivativeCP


First of all, numerical derivatives are quite often needed when dealing with engineering problems. Since derivatives can not be always calculated analytically we need numerical methods.
Furthermore, the CoolProp function returns c_p = 1194.93427991 J/kg, while the own function returns c_p = 1194.93427992 J/kg. This means, the results are identical. Hence, the internal CoolProp function is equal to the own approach. For example, using a forward difference in the own function would lead to a c_p = 1194.93356274 J/kg. This results differs at the third decimal place.


  • Be careful which differencing method you choose close to saturated liquid or saturated vapor line.
  • In the two-phase region c_p tends to infinity. This is because a pure component is not changing its temperature at constant pressure. A temperature glide can occur for mixtures.

What did you think?

I’d like to hear what you think about this post.

Let me know by leaving a comment below.